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section. 
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SECTION  I 


INTRODUCTION  AND  SUMMARY 

This  report  presents  the  results  of  basic  research  in  Model  Algorithmi 
Control  (MAC)  with  application  to  flight  control  systems.  MAC  is  a  new 
digital  control  design  approach  which  relies  on  the  following  features: 

(i)  an  internal  model  of  the  system  to  be  controlled 

(ii)  a  reference  trajectory  description  of  desired  closed  loop 
behavior 

and  (iii)  an  on-line  optimization  of  future  control  inputs  to  produce 
the  desired  performance. 

MAC  is  well  suited  to  the  new  generation  of  microprocessors,  and  is  appli¬ 
cable  to  a  wide  variety  of  aerospace  problems. 

The  purpose  of  this  report  is  to  explore  the  theoretical  foundations 
and  robustness  of  MAC  and  to  examine  its  application  to  flight  control 
problems.  A  missile  attitude  control  simulation  was  chosen  as  a  typical 
control  design  problem  and  used  to  demonstrate  the  application  of  MAC 
theory.  Section  II  presents  a  description  of  the  MAC  design  philosophy 
for  general  systems  and  a  discussion  of  its  application  to  linear  multi- 
variable  problems  in  particular.  Section  III  presents  a  mathematical  for¬ 
mulation  of  an  idealized  discrete-time  MAC  and  examines  robustness,  noise 
behavior,  and  the  effects  of  input  constraints.  Section  IV  extends  the 
theoretical  analysis  of  MAC  to  non-minimum  phase  and  time-delay  systems. 
Section  V  examines  a  continuous-time  MAC  in  order  to  provide  a  more  general 
mathematical  framework  and  application  guidelines  for  the  techniques. 
Section  VI  then  discusses  MAC  design  for  continuous  systems  with  discrete 
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noisy  observations  and  unobserved  outputs,  and  provides  a  formulation  of 
the  problem  of  sample  rate  selection.  Section  VII  discusses  improvements 
to  the  MAC  predictor  to  handle  process  and  measurement  noise. 

Next,  the  application  of  MAC  and  its  software  package  IDCOM  (for 
Identification  and  Command)  to  a  missile  control  example  is  discussed  in 
Section  VIII.  IDCOM  is  described  functionally,  and  the  missile  simulation 
details  are  provided.  The  simulation  results  for  several  different  test 
series  are  given  in  Section  IX.  Section  X  presents  overall  conclusions 
and  recommendations  for  future  work  in  this  area. 

The  major  conclusion  of  this  report  is  that  MAC  is  a  very  flexible, 
powerful  and  promising  technique  for  linear  multivariable  control  design. 
Its  theoretical  properties  of  optimality,  convergence,  robustness  and 
stability  are  verified  by  simulations  of  the  missile  attitude  control 
example.  It  is  recommended  that  future  studies  of  MAC  concentrate  on  its 
extensions  to  the  adaptive  case  and  on  flight  test  demonstration  of  MAC 
concepts. 
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SECTION  II 


MAC  DESCRIPTION 


Overview 

The  purpose  of  this  section  is  to  provide  a  general  introduction  to 
MAC  philosophy  and  to  compare  it  qualitatively  with  other  control  design 
techniques.  The  three  main  principles  of  MAC  are  discussed  in  Section  2.1; 
the  special  features  of  MAC  linear  multivariable  systems  and  a  comparison 
with  other  modern  control  techniques  are  given  in  Sections  2.2  and  2.3. 

The  remaining  two  sections  (2.4  and  2.5)  are  devoted  to  a  component 
functional  description  of  stochastic  MAC  and  a  comparison  with  Self- 
Tuning  Regulators  and  Model  Reference  Adaptive  Control  techniques. 

2. 1  MAC  General  Description 

Some  of  the  recurring  problems  in  the  application  of  modern  control 
theory  to  flight  control  systems  are:  (i)  model  selection,  (ii)  incor¬ 
poration  of  state  and  control  constraints,  and  (iii)  robustness  and  sen¬ 
sitivity  to  unknown  parameters  and  disturbances. 

Here  we  present  a  technique  called  MAC  which  treats  the  above  prob¬ 
lems  in  an  effective  manner.  The  technique  was  originally  developed 
for  industrial  applications  in  France.  It  is  based  on  an  identification- 
optimization  approach,  which  is  very  general  in  nature,  and  is  a  digi¬ 
tal  technique  which  makes  full  use  of  the  capabilities  and  growth 
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potential  of  current  microprocessors.  The  MAC  strategy  relies  on  the 
following  three  features:  1)  internal  model  of  the  system,  2)  reference 
traiectorv  and  outout  constraints,  and  3)  control  traiectorv  computa¬ 
tion. 

2.1.1  Internal  Model  of  the  System 

The  multivariable  system  to  be  controlled  is  represented  by  a  mathe¬ 
matical  model  in  time-domain  of  the  input-output  type  (see  Figure  1). 

For  linear  systems,  the  model  is  of  the  impulse  response  type,  a  repre¬ 
sentation  which  has  certain  distinct  advantages  over  the  state  space 
representation  or  the  transfer  function  representation  for  multivariable 
control.  For  nonlinear  systems,  both  the  state  space  and  tne  input- 
output  representations  have  certain  advantages  and  disadvantages.  In 
applications,  one  may  use  either  or  both  depending  upon  the  nature  of  non- 
linearities  and  the  complexity  of  the  resulting  controller.  The  purpose 
of  the  internal  model  is  to  have  a  flexible  representation  of  the  controlled 
system  stored  in  the  computer  memory,  which  can  be  updated  as  the  system 
changes  and  which  can  be  used  at  any  instant  to  predict  the  future  be¬ 
havior  of  the  system  under  different  control  inputs.  The  internal  model 
of  the  system  is  used  to  compute  optimal  inputs,  to  detect  process  changes, 
sensor  malfunctions  and  severe  faults.  The  inputs  and  current  output 
of  the  internal  model  are  updated  according  to  the  actual  observed  values 
of  these  variables,  but  any  large  difference  between  the  computed  and  the 
actual  values  gives  important  clues  as  to  the  malfunctioning  of  sensors 


and  actuators. 


Figure  1.  Impulse  response  representation  of  a  linear  system. 

The  internal  model  of  the  system  is  generally  obtained  via  off-line 
identification,  using  either  a  physical  model  structure  when  one  is  easily 
available  or  an  input-output  representation  model  such  as  the  impulse  re¬ 
sponse  model,  which  may  change  with  the  operating  point.  Some  form  of 
on-line  parameter  identification  may  also  be  done  in  those  cases  where 
large  random  variations  of  system  parameters  are  expected.  It  has  been 
found  from  experience  that  the  robustness  of  MAC  is  sufficient  to  take 
care  of  small  parameter  changes. 

2.1.2  Reference  Trajectory  and  Output  Constraints 

The  desired  response  of  the  closed-loop  system  is  specified  in  the 

form  of  a  reference  trajectory  and  constraints  are  updated  on-line  using 

the  actual  output  of  the  system.  It  is  possible  to  handle  output-dependent 

constraints^  in  this  fashion  and  to  eliminate  the  steady  state  measurement 

^General  state-dependent  constraints  can  also  be  handled  if  an  observer  (or 
filter)  is  used  to  recreate  the  states  from  the  measured  outputs. 
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offset  error.  It  should  be  noticed  that  the  specification  of  reference 
trajectories  and  constraints  is  much  easier  and  natural  than  the  specifi¬ 
cation  of  a  scalar  performance  index.  Typically,  the  requirements  on 
controlled  outputs  are  stated  as  "no  overshoot,"  "fast  response,"  or 
"within  maximum  and  minimum  limits."  These  requirements  are  difficult  to 
express  in  a  scalar  quadratic  performance  index,  but  they  are  easily  con¬ 
verted  into  desired  trajectories  and  constraints  which  can  also  be  ex¬ 
plained  to  the  operators  without  any  difficulty.  In  a  hierarchtcal  control 
system,  reference  trajectories  for  lower  levels  may  be  specified  by  the 
higher  levels  through  an  optimization  process,  since  the  controlled  outputs 
of  lower  levels  are  the  control  variables  for  the  higher  levels.  The  ref¬ 
erence  trajectories  at  the  highest  level  would  generally  be  obtained  via 
a  combination  of  on-line  and  off-line  optimization. 

The  concept  of  using  reference  trajectories  is  more  general  than 
model -fol lowing.  Firstly,  it  may  not  be  possible  to  represent  a  reference 
trajectory  by  a  simple  model,  and  secondly,  under  sensor  or  actuator 
fault  conditions,  one  may  have  to  relax  the  system  requirements  to  con¬ 
trol  within  a  band  or  tolerance  limit.  Certain  control  problems  in¬ 
volving  more  controlled  outputs  than  control  variables  are  formulated 
more  correctly  as  band-control  problems  rather  than  model-following  or 
scalar  performance  index  problems. 

2.1.3  Control  Trajectory  Computation 

Controls  are  computed,  in  general,  for  a  number  of  future  time  points 
using  an  iterative  optimization  technique  which  minimizes  the  distance 
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between  the  desired  reference  trajectory  and  output  trajectory  predicted 
by  the  internal  model,  while  keeping  all  the  output,  state  and  control 
constraints  satisfied.  The  complexity  of  the  control  algorithm  is  directly 
dependent  on  the  structure  of  the  internal  model,  the  number  of  inputs 
and  outputs  and  the  constraints.  For  linear  systems,  the  impulse  response 
representation  results  in  a  simple,  fast,  projection-type  quadratic  pro¬ 
gramming  solution  which  can  be  implemented  in  micro/mini -computers  of 
the  present  generation.  The  actual  dimensionality  of  the  state  does  not 
increase  the  complexity  of  the  algorithm  as  it  would  in  a  state  vector 
representation. 

2.2  Special  Features  of  MAC  and  IDCOM  (IDentification  and  COMmand)  for 

Linear  Multivariable  Control 

The  MAC  design  approach  has  been  found  to  be  flexible  and  robust. 

It  is  well  suited  for  the  evolving  microprocessor  technology  providing 
high  speed  memory  and  fast  computation  times  for  basic  calculations  such 
as  convolutions.  The  universality  of  the  impulse  response  representation 
leads  to  a  unified  design  approach  for  systems  of  all  orders.  Further¬ 
more,  the  parameter-linearity  of  this  representation  leads  to  a  duality 
between  identification  and  control. 

MAC  is  implemented  by  a  computer  program  called  IDCOM.  The  special 
features  of  IDCOM  are:  (i)  no  model  order  reduction  is  required  since  an 
inpulse  response  representation  is  used,  (ii)  Input  magnitude  and  rate 
constraints  are  handled  directly  and  exactly,  (iii)  The  control  law  is 
time-varying  and  the  closed  loop  response  is  robust  to  parameter  changes, 
(iv)  Gain  scheduling  is  replaced  by  on-line  updating  of  the  internal  model 
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using  operating  data  and  parameter  estimation  techniques,  thereby  reducing 
reliance  on  theoretical  models  of  the  system,  (v)  The  same  algorithm  is 
used  for  impulse  response  identification  and  for  control  law  computation, 
thereby  simplifying  the  hardware  requirements,  (vi)  The  control  laws 
can  be  modified  on-line  in  case  of  sensor  failures  or  degraded  system 
performance. 

2.3  Comparison  of  IDCOM  with  Other  Modern  Control  Techniques 

IDCOM,  like  most  other  modern  control  techniques,  is  a  multivariable 
time-domain  technique.  It  is  different,  however,  in  several  important 
ways.  These  differences  are  important  since  they  account  for  the  success 
of  IDCOM  in  control  applications  where  the  lack  of  good  models  has  pre¬ 
vented  the  use  of  other  modern  control  techniques  (Richalet  et  al^.,  1978). 
The  important  differences  are: 

(i)  Robustness.  Modern  control  techniques  generally  use  full  state 
feedback  based  on  a  parametric  state  vector  model.  IDCOM,  on  the  other 
hand,  uses  output  feedback  based  on  an  impulse  response  model.  It  can 
be  argued  that  the  effect  of  modeling  errors  on  IDCOM  will  be  less  com¬ 
pared  to  that  on  state  vector  techniques. 

(ii)  Implementation.  Modern  control  techniques  based  on  state  vector 
models  generally  require  the  solution  of  matrix  Riccati  equations.  Since 
these  equations  are  computationally  time-consuming  to  solve  for  systems  of 
high  order,  the  practical  implementation  of  these  techniques  requires 
model  reduction  and  off-line  solution  of  the  Riccati  equation.  In  prac¬ 
tice,  only  the  steady  state  gains  from  the  Riccati  equation  are  imple¬ 
mented  and  these  gains  are  scheduled  as  a  function  of  the  operating  point 
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(e.g.,  flight  condition).  This  implementation  does  not  allow  for 
on-line  changes  in  system  model  and  performance  criteria.  The  IDCOM 
approach,  on  the  other  hand,  is  much  more  flexible  and  adaptive,  since 
the  model,  criteria,  and  sampling  rates  can  be  adjusted  on-line.  This 
flexibility  comes  from  the  use  of  the  impulse  response  representation, 
which  has  the  additional  advantages  of  conceptual  simplicity,  ease  of 
identification  and  elimination  of  the  model  order  reduction  problem. 

(iii)  Criterion  Function  and  Constraints.  Most  of  the  modern 
control  techniques  require  specification  of  a  sca’ar  criterion  function. 

In  practice,  this  is  very  difficult  or  impossible  to  do.  The  specifica¬ 
tion  of  the  information  used  in  IDCOM  is  easier  and  more  natural.  Basi¬ 
cally,  the  specification  of  desired  responses  and  constraints  is  simpler 
than  the  specification  of  a  scalar  performance  index.  It  is  also  possible 
to  set  up  a  priority  list  on  the  control  of  outputs  in  IDCOM  and  even  to 
make  them  conditional  on  the  occurrence  of  some  future  events.  For  example, 
the  failure  of  a  sensor  can  be  detected  and  the  control  strategy  can  be 
shifted  accordingly.  The  problems  of  sensor  blackout  or  computer  over¬ 
load  can  also  be  handled  since  controls  are  computed  and  stored  for 
several  future  time  points  and  these  can  be  used  till  new  information  or 

computational  capability  becomes  available.  The  exact  satisfaction  of 
control  and  output  constraints  is  absolutely  essential  in  many  applications. 
Such  constraints  are  handled  much  more  easily  in  IDCOM,  than'  s  to  the 
impulse  response  representation  of  the  system. 

(iv)  Duality  of  Identification  and  Control.  From  a  practical  stand¬ 
point,  the  duality  of  identification  and  control  for  an  impulse  response 
model  is  much  more  useful  than  the  duality  between  estimation  and  control 


for  state  vector  models.  The  former  can  be  used  directly  to  increase  the 
robustness  of  the  controller  by  on-line  identification  whenever  the  resi¬ 
duals  between  the  model  output  and  the  system  output  are  consistently 
outside  some  prespecified  confidence  bands. 

Some  of  the  modern  control  approaches  which  are  closer  to  MAC  are  the 
Model  Reference  Adaptive  (MRA)  (Landau,  1974)  and  the  Self-Tuning  Regula-or 
(STR)  (Astrom,  1980)  approaches.  They  are,  however,  not  exactly  similar 
since  the  specification  of  reference  models  and  the  computations  of  con¬ 
trol  are  done  differently,  as  discussed  in  Section  2.5  below. 

2.4  Component  Functional  Description 

In  order  to  explain  the  operation  of  MAC  and  compare  it  to  other 
approaches,  it  is  useful  to  describe  the  functions  of  each  explicit  MAC 
component.  Figure  2  shows  different  elements  of  a  fairly  general  MAC 
scheme  consisting  of  seven  different  functional  blocks.  These  are: 

1.  Internal  model  of  the  plant  which  is  meant  to  be  a  close  repre¬ 
sentation  of  the  system  (this  model  may  be  in  state  vector  form,  linear 
or  nonlinear,  impulse  response  form,  transfer  function  form  or  even  table 
look-up,  depending  on  the  application). 

2.  Reference  trajectory  generator,  which  uses  the  latest  plant  out¬ 
puts  and  desired  set  point  conditions  to  generate  desired  reference  tra¬ 
jectories  for  the  controlled  closed  loop  system  to  follow. 

3.  Control  signal  generator  performs  most  of  the  main  computations 
within  MAC  to  produce  control  signals  for  the  actuators.  The  inputs  to 
this  block  consist  of  the  predicted  error  between  the  output  of  the 
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Figure  2.  Relation  of  adaptive  MAC  to  other  adaptive  © 
control  schemes .  (% 
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reference  trajectory  generator,  ©,  and  the  output  predictor,  ®,  along  with 
the  control  and  state  constraints, 

The  control  algorithm  is  iterative  involving  several  trials  of  dif¬ 
ferent  control  inputs  to  calculate  the  best  input  for  minimizing  the 
tracking  error  without  overstressing  the  actuators  and  the  computational 
facilities  of  the  process  control  computer.  The  model  used  for  control 
computations  may  not  be  the  same  as  the  internal  model ,  (I),  of  the  plant. 

4.  Identifier/ scheduler  performs  the  function  of  updating  the  inter¬ 
nal  model  0  of  the  plant  under  changing  operating  and  plant  conditions. 

If  the  changes  in  the  plant  can  be  expressed  as  a  deterministic  function 
of  the  operating  state,  off-line  or  periodic  identification  followed  by 
scheduling  will  be  adequate. 

However,  if  the  changes  in  the  plant  are  rapid  and  random,  on-line 
identification  is  required,  which  may  be  passive  or  active.  In  passive 
on-line  identification,  no  extra  signals  are  injected  into  the  system  to 
enhance  identifiabil ity  or  to  improve  the  speed  of  identification.  In 
active  on-line  identification,  extra  input  signals  are  used  to  enhanc 
identifiabil ity  in  order  to  match  the  frequency  of  adaptation  with  the 
frequency  of  plant  variations  (see  Mehra  et_  al_. ,  1978). 

Input  signal  design  techniques  (Mehra,  1980)  play  an  important  role 
in  active  adaptation.  The  identifier  can  also  be  augmented  to  perform 
Fault  Detection,  Diagnosis  and  Prognosis  (DDP)  functions. 

5.  Control,  state,  and  measurement  constraints,  which  are  generally 

of  the  amplitude  and  rate  type,  may  in  some  cases  (e.g.,  mixed  state  and 

control  constrain  s)  be  much  more  complicated. 

2 

Circled  numbers  refer  to  blocks  in  Figure  2. 
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6.  Output  predictor  is  an  important  element  of  MAC  since  the  accuracy 
of  plant  predictions  contributes  directly  to  the  performance  of  any 
control  system.  The  predictor  may  be  of  the  state  type  (Kalman,  extended 
Kalman,  etc.),  Luenberger  observer  type,  impulse  response,  or  ARMA  (Box 
and  Jenkins,  1976  or  Astrom,  1970).  The  latter  type  is  common  in  Self- 
Tuning  Regulator  (STR)  and  Model  Reference  Adaptive  (MRA)  control 
applications.  They  are  also  discussed  in  Section  VII. 

7.  Deconvolution  generator  is  used  in  those  cases  where  a  direct 
use  of  the  internal  model  will  lead  to  undesirable  control  behavior, 
e.g.,  nonminimum  phase  or  time  delay  systems.  Several  techniques  for 
generating  deconvolution  models  are  discussed  in  Section  IV  and  Appendix  A. 

2 . 5  Comparison  with  STR  and  MRA 

In  Figure  2,  a  comparison  of  MAC,  STR,  and  MRA  is  given  in  terms  of 
function  blocks  used  in  each  scheme.  For  example,  in  STR  with  implicit 
plant  model,  functions  OD,  3,  and  i4)  are  combined  into  one  block,  i.e., 
the  controller  directly  identifies  the  gains  of  the  control  law  without 
explicit  identification  of  the  plant  model. 

In  this  scheme,  blocks  3,  '51,  §)  and  X  are  either  absent  or  they  are 
combined  with>;T,  3  and  Such  a  combination  results  in  a  simple  design 
at  the  expense  of  flexibility.  For  example,  one  cannot  change  on-line 
reference  trajectories  or  account  for  failed  sensors.  In  explicit  STR, 
an  ARMA  model  is  identified  on-line  and  then  used  to  develop  the  control 
law.  This  scheme  is  more  suitable  for  nonminimum  phase  systems,  but  still 
cannot  handle  constraints  exactly. 


MRA  may  be  implicit  or  explicit  with  respect  to  the  reference  model. 
The  difference  is  indicated  in  Figure  2,  where  ®  +  ®  implies  that  func¬ 
tions  ©  and  ©  are  combined  in  the  implicit  MRA  technique. 

The  obvious  difference  between  MAC  and  the  other  techniques  is 
separation  of  different  functions.  In  some  practical  applications,  all 
the  functions  may  not  be  required,  but  the  great  generality,  robustness 
and  flexibility  of  MAC  come  from  having  all  these  functions  in  the  soft¬ 
ware.  With  modern  advances  in  digital  hardware,  one  incurs  very  little 
cost  for  having  separate  functions,  whereas  the  payoff  in  terms  of  per¬ 
formance  and  flexibility  of  the  control  scheme  can  be  quite  high. 

2.6  Summary 

In  this  section,  the  philosophy  and  basic  concepts  of  MAC  were  dis¬ 
cussed.  The  similarities  and  differences  between  MAC  and  other  modern 
control  techniques  were  outlined,  and  the  generality  and  flexibility  of 
MAC  for  control  design  purposes  were  emphasized. 
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SECTION  III 


MATHEMATICAL  FORMULATION  AND  ROBUSTNESS  OF  IDEALIZED  DISCRETE-TIME  MAC 

Overview 

The  puroose  of  this  and  some  of  the  subsequent  sections  is  to  develop  a 
mathematical  framework  for  an  analytical  study  of  MAC.  The  description  of 
MAC  in  Section  II  reveals  that  a  great  deal  of  intuition  is  involved  in  MAC 

3 

design.  Therefore,  a  complete  mathematical  analysis  of  MAC  properties  is 
extremely  difficult.  Several  simplifying  assumptions  are  required  to  make 
the  analysis  tractable.  The  purpose  of  mathematical  analysis  is  to  serve  as 
a  guide  for  MAC  design  rather  than  a  replacement  for  it.  In  this  section, 
we  analyze  unconstrained  linear  single-input  single-output  discrete-time, 
minimum  phase  systems.  Both  deterministic  and  stochastic  (colored  output 
additive  noise)  systems  are  considered.  The  robustness  problem  is  discussed 
and  performance  "measures"  of  robustness  are  proposed.  Some  of  results 
presented  here  are  more  generally  applicable  to  optimizing  type  of  control 
laws  of  which  MAC  is  one  example.  In  particular,  a  close  link  is  estab¬ 
lished  between  the  stability  and  closed  loop  properties  of  optimizing  type 
and  feedback  type  of  control  laws. 

The  organization  of  this  section  is  as  follows.  Section  3.1  presents 

an  overview  of  the  major  functions  of  MAC  (prediction-optimization). 

Section  3.2  treats  the  case  of  a  single  input-single  output  system  where 

the  inputs  are  free  of  constraints,  and  introduces  linear  closed  loop  and 

open  loop  predictors.  Section  3.3  is  concerned  with  the  stability  analysis, 
3  ......  . 

The  reader  may  be  interested  in  looking  at  Sections  8.1,  8.2  and  8.3  for 
further  details  of  the  MAC  algorithm. 
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both  from  a  time  domain  and  frequency  domain  viewpoint.  Section  3.4 
investigates  the  effect  of  an  output  additive  noise  and  compares  the  per¬ 
formances  of  the  control  with  respect  to  reference  trajectories  under 
open  and  closed  loop  prediction  schemes.  Section  3.5  poses  the  robust¬ 
ness  problem  and  its  relation  to  stability.  Section  3.6  defines  two  per¬ 
formance  indices  for  robustness  analysis.  Finally,  in  Section  3.7  con¬ 
straints  on  input  are  introduced  and  the  analysis  of  the  constrained  system 
is  related  to  the  analysis  of  the  unconstrained  one. 

3.1  Mathematical  Formulation 

As  discussed  in  Section  2.5,  MAC  is  conceptually  similar  to  a  model 
reference  adaptive  type  of  control  with  some  important  differences  in 
practical  implementation.  It  involves  (i)  dynamic  models  for  system  rep¬ 
resentation  and  prediction,  (ii)  a  reference  trajectory  and  (iii)  an 
optimality  criterion  leading  to  the  optimal  control .  A  schematic  rep¬ 
resentation  of  MAC  (Figure  3)  displays  the  main  operations  performed: 

(i)  prediction  of  the  future  output  for  T  steps  ahead  based  on  the  output 
(y ( t ) )  of  the  actual  plant  and  on  the  computed  output  (yM(t))  of  the 
plant's  model,  (ii)  calculation  of  the  future  reference  trajectory  based 
on  the  actual  output  (y ( t ) )  of  the  plant  and  a  desired  set  point  c(t), 
and  (iii)  computation  of  the  input  u(t)  based  on  the  above  predictions  and 
using  a  certain  optimization  criterion. 

3.1.1  Representation  and  Prediction 

The  system  is  represented  by  its  impulse  responses,  the  identification 
of  which  can  be  done  both  on-line  and  off-line.  However,  in  most  cases 
the  off-line  identification  is  accurate  enough  for  the  purpose  of  control 
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Figure  3.  Model  Algorithmic  Control 


and  one  can  avoid  the  cost  and  complexity  of  an  on-line  identification 
procedure.  This  is  due  to  the  particular  redundancy  of  the  impulse  response 
representation  which  allows  a  considerable  enhancement  of  the  robustness  of 
the  control  scheme  against  identification  errors  and  parameter  pertur¬ 
bations  (Mehra  et  aK  ,  1979;  Mereau  et  al . ,  1978). 

The  formal  representation  of  the  system  is  as  follows: 

y(t+l)  =  hTu(t)  =  h0u(t)  +  h1u(t-l)  +  ...  +  hNu(t-N)  (3.1) 

where  y(t+l)  is  the  plant  output  at  time  t+1;  hJeIRN+1  denotes  the  plant 
impulse  response;  u(t-j)eft>-]R  for  j  =  0,...,N.  u(t-j)  is  the  input  at 
time  (t-j)  to  the  plant,  ft  is  the  constraint  set  of  the  input. 

The  model  (3.1)  is  also  called  the  Actual  Model  to  emphasize  that  hT 
represents  the  actual  impulse  response.  But  since  such  perfect  knowledge 
V)}£  denotes  the  sequence  {c( t+1) , . . .  ,c(t+T) } . 
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of  the  plant  impulse  response  is  usually  not  possible,  one  has  to  use  an 
approximation  to  h^.  The  model  corresponding  to  this  latter  impulse 
response  is  then: 

yM(t+i)  =  hTu(t)  (3.2) 

The  above  model,  together  with  the  past  history  of  the  plant  output 
denoted  by  Y(t)=  {y(-r),  t s  t}  is  used  to  predict  the  future  value  of  the 
output.  Various  prediction  schemes  are  conceivable.  In  this  section  we 
will  limit  ourselves  to  simple  open  loop  and  closed  loop  prediction. 
Figure  4,  where  yp(t)  refers  to  the  predicted  value  of  the  output, 
summarizes  the  representation-prediction  part  of  MAC. 


Figure  4.  MAC  prediction. 


3.1.2  The  Reference  Trajectory 

The  purpose  of  the  control  is  to  lead  the  output  y(t)  along  a  desired, 
and  generally  smooth,  path  to  an  ultimate  set  point  C.  Such  a  path  is 
called  a  reference  trajectory.  In  the  present  section  the  reference  tra¬ 
jectory  is  of  first  order  and  is  initiated  on  the  output  of  the  plant  at 
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time  t q : 

yr(tn+k)  =  akyr(t0)  +  (l-ak)c  k=l,...,T 
r  °  r  U  (3.3) 

yr(t0)  =  y(tQ) 

The  reference  trajectory  can  be  chosen  to  be  of  higher  order  and  the  set 
point  c  can  be  made  time  varyinq.  These  latter  extensions  would  not 
bring  any  conceptual  difficulty  in  the  solution  of  the  control  problem, 
although  the  actual  computation  of  the  inputs  becomes  more  complex,  but 
not  intractable,  as  the  order  of  the  reference  trajectory  increases. 

3.1.3  The  Optimality  Criterion  and  the  Optimum  Control  Strategy 
The  optimality  criterion  should  reflect  the  previously  mentioned 
purpose  of  following  the  reference  path  to  the  desired  set  point  c. 

That  can  be  done  by  defining  the  optimum  control  strategy  as  the  one 
which  minimizes  over  a  certain  horizon  in  the  future,  the  deviation  of 
the  predicted  outputs  from  the  reference  path.  Fornelly,  at  each  instant 
tQ,  the  optimum  set  of  T  future  inputs  {u *(tQ) ,  u *(tQ+l) , . . . ,  u  *(tQ+T-l)} 
are  such  that  the  predicted  T  outputs  {yp( tQ+l) ,. . . ,yp(t0+T)}  are  as 
close  as  possible,  in  the  sense  of  a  weighted  Euclidean  norm,  to  the  ref¬ 
erence  trajectory  yr-  Therefore  the  function  to  minimize  is: 

JT  =  I  (yp(t+k)  -  y  (t+k))Z  w.  (3.4) 

1  k=l  K  r  K 

w^  :  nonnegative  weighting  factor 

Note  that  at  time  t^,  the  determination  of  T  optimal  inputs  u*(tg+k) 
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(k  =  0,. . .  ,T-1)  is  done  by  solving  a  static  optimization  problem: 


Minimize  J-j.(u(tg) , . . .  ,u  (tQ+T-l) ) 
s.t.  u  (tQ+k)  e  Q.  for  k=  0,. ..  ,T-1 


(3.5) 


Various  types  of  algorithms  can  be  used  to  solve  (3.5)  and  to  determine  the 
set  U*(tQ)  =  {u  *(tg) , . . . ,u*(tQ+T-l ) }  (Luenberger,  1  973).  The  interval  {tQ,tg+T} 
is  called  the  horizon  of  control  evaluation,  and  sometimes  the  horizon 
of  prediction  since  at  time  tg  one  has  to  predict  a  set  of  T  outputs 
yp(t+k),  (k=l»...,T).  Once  the  set  U*(tg)  of  the  T  optimum  inputs  is 
determined,  it  is  possible  to  wait  up  to  T  periods  before  observing  the 
actual  output  y(t),  reinitializing  yr(t),  predicting  yp‘s  and  computing 
the  next  set  U*(tg+T).  This  means  that  all  the  elements  of  the  optimum 
control  set  U*(tQ)  have  been  actually  applied  to  the  plant. 

A  more  appealing  strategy  consists  of  applying  the  first  few  optimal 
inputs  u*(t+k),  (k  =  0,l,...,p  with  p«T)  before  reinitializing  y^  and 
computing  the  next  T  optimum  inputs.  In  the  limit,  if  the  computing 
facilities  allow,  one  would  apply  only  the  first  optimum  input  u*(t)  of 
the  set  U*(t),  observe  y(t+l),  initialize  the  reference  trajectory  yr  on 
the  observed  value  y(t+l)  and  solve  the  optimization  problem  (3.5)  at 
time  t+1.  In  this  latter  case  the  unused  optimum  inputs  u*(t+l),..., 
u*(t+T-l)  of  U*(t)  (computed  at  time  t)  will  serve  as  starting  values 
for  the  numerical  algorithm  determining  U*(t+1),  that  is,  solving  (3.5) 
at  time  (t+1).  Only  this  latter  case  will  be  considered  in  this  section. 

Figure  5  displays  the  optimum  inputs  and  the  corresponding  predicted 
output  over  the  horizon  T.  To  visualize  the  overall  control  procedure 
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corresponding  to  the  last  mentioned  case,  one  must  repeat  the  same  picture 

for  t+1,  t+2, -  Fiqure  6  is  a  block  reDresentation  of  the  whole  control 

scheme. 


Figure  5.  .'1AC  inputs  and  outputs. 


* 


Figure 


6. 


Control  block  representation. 


3.2  Unconstrained  Input 

>  . 

3.2.1  Linear  Prediction 

So  far  the  characterization  of  the  prediction  scheme  has  been  kept  general. 
Let  us  constrain  the  class  of  predictors  to  be  linear.  Notwithstanding  the 
linearity  of  the  plant,  the  above  class  of  predictors  is  a  practically  con¬ 
venient  choice.  Moreover,  it  simplifies  substantially  the  determination 
of  the  optimum  control  sequence.  Indeed,  the  linear  character  of  yp(t) 
causes  the  function  Jy  to  be  quadratic  in  the  inputs  u(t),  and  hence  to 
become  a  candidate  for  fast  quadratic  optimization  algorithms. 

3.2.2  One  Step  Prediction  (T  =  1) 

The  assumption  that  the  input  vector  u(t)  is  free  of  constraints, 

M 

i.e.,  n=R  ,  results  in  a  significant  simplification  in  the  optimum  control 
determination.  That  is,  the  length  T  of  the  horizon  of  prediction  does 
not  affect  the  optimum  value  of  the  first  input  u*(t)  to  be  applied.  In 
other  terms,  the  minimization  of  Jy  and  Jy,  ,  with  T1  t  T,  will  result  in 
the  same  first  input  u*(t)  of  the  optimum  sequences  of  length  T'  and  T. 

In  particular,  the  first  element  of  the  sequence  {u* ( t ) ,. . . ,u*(t+T-l)} 
minimizing  Jy  is  identical  to  the  input  minimizing  Jy  .  This  is  an  im¬ 
portant  simplification,  since  it  reduces  a  T-variable  minimization  at 
each  step  to  a  one-variable  minimization.  This  property,  which  is  based 
essentially  on  the  principle  of  superposition  of  linear  systems,  is 
easily  established,  as  the  proof  of  the  following  proposition  demonstrates. 

Proposition  3.1.  If  the  input  to  the  plant  defined  by  (3.1)  is  free  of 
constraints,  the  minimization  of  Jy  (3.4)  and  the  minimization  of  Jy  at  time 
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tQ,  result  in  the  same  first  input  u*( tg) .  That  is,  the  first  element  of 
the  set  U*(tg)  corresponding  to  the  minimization  of  Jy  is  identical  to  the 
input  minimizing  J1 . 

Proof 

T  2 
JT  =  .1  H Wj)  '  wj’  wj>0 

J  1 

The  value  zero  for  Jy  is  achieved  if  and  only  if: 

yp(t0+j)  =  yr(tQ+j)  j  =  1 , —  ,T  (3.6) 

The  right-hand  side  of  equation  (3.6)  depends  on  the  observed  value 
y(tQ),  a  and  c  through  (3-3),  and  hence  is  known  independently  of  the 
input  u(tQ) ,u(tQ+l) ,. . . .  The  left  hand  side  yp(t0+j)  depends  linearly 
°n  yM(tQ+j)  and  on  the  observed  value  {y(i),  T<tg+j-l}.  Therefore,  the 
inputs  u  (tQ)  ,u  (tQ+l) ,... ,u(tQ+j-l)  enter  linearly  in  the  expression 
of  yp(tQ+j).  The  determination  of  the  optimum  input  sequence  is  done  as 
follows.  For  j  =  l  the  equation 

yp(t0+1)  =  yr(Vi)  (3-7^ 

determines  uniquely  u*(tg).  The  next  equation,  corresponding  to  j  =  2: 
yp(t0+2)  =  yr(t0+2) 

in  which  the  value  of  u(tQ)  is  set  equal  to  the  optimum  u*(tQ),  deter¬ 
mines  uniquely  u*(tg+l).  Hence  the  optimal  sequence  { u*  ( tg) ,u*(tQ+T-l)} 
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is  determined  sequentially.  The  existence  and  uniqueness  of  a  solution 
for  each  equation  of  (3.6)  is  assured  because  of  the  linearity  of  yp(t+j). 
Now  let  us  consider  the  cost  J,: 

x 

Ji =  I!  W1)  2 

Clearly,  the  minimum  of  Jy  is  achieved  for: 

yp(t0+l)  -  yr(tg+l)  (3.8) 

which  is  identical  to  (3.6),  resulting  in  the  same  solution  u*(tg).  ■ 

With  respect  to  the  above  proof  we  should  mention  that: 

(1)  The  linear  character  of  yp(t+j)  is  a  sufficient  condition  for 
the  proof  to  hold.  If  the  predictor  yp(t+j)  is  nonlinear  and  its  depen¬ 
dence  on  the  past  inputs  and  the  observed, data  is  denoted  via  a  nonlinear 
operator  F: 

yp(t+j)  =  F ( u ( t+ j - 1 ) ,  u ( t+ j - 2 ) .... ;  y(i ) ,  r<t+j-l) 

then  the  necessary  general  condition  for  the  proof  of  proposition  3.1  to 
hold  is  that 

F(  • ,  u  ( t+j  -2 ),...,  u  ( t+ j  -N- 1 ) ;  y(x) ,  x<t+j-l) 

must  be  a  mapping  of  R  onto  R.  In  the  asc  of  a  nonlinear  predictor,  the 
optimal  input  sequence  is  generally  not  unique. 

(2)  The  T  optimal  inputs  u*( tg), . . . ,u*(tg+T-l)  obtained  by  minimizing 
Jy  are  not  identical  to  the  input  sequence  obtained  by  T  consecutive  mini¬ 
mizations  of  Jy.  This  is  partly  because  of  the  right-hand  side  of 
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equation  (3.6).  Indeed,  while  minimizing  the  distance  over  a  horizon  of 


length  T,  yr(tQ+j)  is  initialized  on  y( tQ)  for  the  whole  length  of  in¬ 
terval  T;  but  in  the  case  of  T  consecutive  minimizations  of  at  time  tQ, 
tj+1,. . . ,tg+T-l,  the  reference  trajectory  is  continuously  initialized 
on  the  observed  y(tQ) ,y(tQ+l) ,. . . ,y(tQ+T-l);  that  is: 

yr(tQ+j)  =  xy ( tQ+ j - 1 )  +  (l-ot)c  for  all  j=l,...,T 

3.2.3  Open  L oop  and  Closed  Loop  Prediction 

In  this  section  two  simple  one  step  prediction  (T=  1)  schemes,  one 
open  loop  and  the  other  closed  loop,  are  investigated.  The  input  is 
assumed  to  be  free  of  constraints. 

The  simplest  one  step  open  loop  predictor  that  one  can  imagine  is: 

yp(t+l)  =  yM(t+l)  =  hTu(t)  (3.9) 

that  is,  the  model  of  the  plant  is  used  to  nredict  the  outout  one  step 
ahead.  The  main  inconvenience  of  this  open  loop  predictor  is,  as  one 
might  expect,  that  the  output  y(t)  of  the  controlled  plant  would  not 
converge  to  its  desired  final  value  c.  This  is  easily  established  as 
follows.  The  optimum  input  sequence  u(t),  in  the  sense  of  minimizing 
Jp  satisfies  (3.8),  which  together  with  (3.3)  and  (3.9)  implies: 

yM(tH)  -  hTu(t)  =  yp(t+l)  =  nty(t)  +  (l-i)c  (3.10) 


Now  assuming  that  the  above  recursive  equation  is  convergent  (i.e., 
the  system  is  stable),  the  equilibrium  input  u*  and  plant  output  y*  are 
given  as: 


( l-a)c 
hTl  -  cthTl 


a  constant  scalar 


y*  = 


hTi 


-U-a)c 


hTl  -  ahTl 


a  constant  scalar 


1  =  (1,1,...  ,1)  e  IR 


N+l 


(3.12) 


Since  in  general  the  model  impulse  response  is  different  from  the 
plant  one  hT,  the  ultimate  value  y*  differs  from  the  desired  value  c.  It 
is  seen  that  the  necessary  and  sufficient  condition  resulting  in  a  nonbiased 
open  loop  prediction,  i.e.,  y*~c,  is  hTl=hTl.  But  such  a  condition  is 
more  of  academic  interest  and  has  little  chance  of  holding  in  practical 
situations.  However,  it  is  important  to  note  that  the  bias 


'T  T  c 
h  i  -  ah  1 


can  be  estimated  since  h^  is  known  and  h^l,  i.e.,  the  gain  of  the  process, 
can  be  estimated  from  the  step  response  of  the  plant. 

To  overcome  the  above  bias  problem,  one  uses  a  closed  looo  orediction, 
a  sinnle  case  of  which  is: 


yPU+D  *  yM(t+l)  +  (y(t)-yM(t)) 
=  y(t)  +  (yM(t+l)  - yM( t) ) 


(3.13) 


where  (y(t)-yM(t))  represents  a  correction  term  which  is  similar  to  the 
innovation  term  of  a  Kalman  filter.  Then  the  optimum  control,  by  (3.7),  is 


such  that: 


ay(t)  +  (l-a)c  =  y(t)  +(yM(t+l)  -  yM( t) ) 

Now  letting  t -*•“>,  and  assuming  that  the  system  is  stable  (i.e., 
one  deduces: 

lim  y(t)  =  c 

t-+oo 


(3.14) 


lim  y(t)  <  °°), 

t*+oo 


that  is,  the  desired  set  point  c  is  reached  asymptotically.  The  previous 
bias  has  vanished. 


3.3  Stability 

In  the  previous  section,  it  has  been  assumed  that  the  optimally 
controlled  system  (both  with  open  loop  and  closed  loop  prediction)  is 
stable.  This  section  addresses  the  stability  condition  both  in  time  do¬ 
main  and  in  z-transform  domain. 

In  both  cases,  the  optimum  sequence  of  inputs  u*(t)  ,is  generated  by 
an  autoregressive  equation,  which  results  from  (3.8)  where  yp(t+l)  and 
y  (t+1)  are  expressed  in  terms  of  inputs: 

open  loop:  (from  equations  (3.9)  and  (3.10)) 

N  N-l 

hnu*( t)  =  a  l  h.u*(t-j-l)  -  l  fi,u*(t-j)  +  (l-a)c 
U  j=0  J  j=l  J 

y(t)  =  hTu*(t-l) 
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(3.15) 


closed  loop:  (from  equations  (3.9)  and  (3.14)) 

N-l  __  N-l  _  N-l 

hnu*(t)  =  l  h,u*(t-l-j)  -  l  h.u*(t-j)  +  (l-a)[c  -  l  h.u*(t-j)] 

U  j=0  J  j=l  3  j=0  J 

(3.16) 

y(t)  =  hTu*(t-l) 

Obviously,  if  the  sequence  u*(t)  tends  to  an  equilibrium  value  u*,  that 
is,  if  the  corresponding  autoregressive  model  is  stable,  then  the  output 
y(t)  tends  to  an  equilibrium  value  y(°°)  which  equals  c  in  the  case  of 
closed  loop  prediction  but  differs  from  c  for  open  loop  prediction.  How¬ 
ever,  the  converse  is  not  true;  that  is,  theoretically  one  may  have  a 
converging  output  y(t),  while  the  input  u*(t)  diverges.  Intuitively  it 
is  clear  that  even  though  |u*(t)|  might  increase  indefinitely,  the  linear 
function  hTu*(t)  may  remain  finite  (at  least  in  theory).  This  is  basically 
what  happens  in  nonminimum  phase  systems.  Such  behavior  is  not  acceptable 
in  applications  because  of  the  nonrealizability  of  infinite  inputs  (Astrom, 

1970;  Astrom  and  Wittenmark,  1974). 

The  boundedness  condition  for  the  input  sequence  u*(t)  is  identica’ 
to  the  stability  of  the  autoregressive  models,  that  is,  the  polynomials: 

z[hQ  +  RjZ-1  +  ...  +  hN_1z_N_1Il  -  a[hQ  +  hjZ  1  +  .. .  +  h^z  N  1]  =  0 

(open  loop)  (3.17) 

and 

(z-l)[h0  +  h1z-1+...+hN_1z'N-1]  +  (l-a)[h0+h1z-1+...+hN_1z'N-1]  =  0 

(closed  loop)  (3.18) 


must  have  all  their  roots  within  the  unit  circle. 
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So  far  we  have  mainly  focused  on  the  relationship  between  the  output 
y(t)  and  the  optimum  input  sequence  u*(t).  Now  let  us  turn  our  attention 
to  the  response  of  an  optimally  controlled  system  to  the  set  point  c,  that 
is,  the  relationship  between  y(t)  and  c(t).  In  terms  of  the  z-transforms , 
one  deduces  from  (3.11)  and  (3.14)  with  some  manipulations: 


open  loop: 


_ (l-jO_ _ 

H(z)  -  :iz  ^H(z) 


(3.19) 


=  z  JH(z)(l-~.)  __ 

H(Z  )  -  uZ  hl(z) 


(3.20) 


closed  loop: 


(l-o 


z'1(l-a)H(z)  +  (1-Z-1)H(Z) 

_ Z~1H(z)(l-Q _ 

Z_1(l-u)H(z)  +  ( l-Z_1)H(z) 


(3.21) 


(3.22) 


with  u(*)»  Y(*)<  H(*)«  H(*)  and  C(*)  denoting  respectively  the  z-transforms 
of  u(t),  y(t),  h(t),  h(t)  and  c(t).  Figures  7  and  8  display  the  open- 
loop  and  the  closed-loop  cases. 


Figure  7.  Open  loop  prediction. 
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Figure  8.  Closed  loop  prediction. 


Let  us  first  note  that  for  a  perfect  identification  of  h,  that  is, 
for  H(z)  =  H(z),  the  open  loop  and  closed  loop  transfer  functions  become 
identical  and  given  by 


uU)  =  (l-g) 

^  (l-az_1)H(z) 

Y(z)  =  1j2 

C(2)  Z-a 


(3.23) 

(3.24) 


Note  that  under  such  perfect  identification  the  transfer  function  of  Y(z) 
with  respect  to  C(z)  is  of  first  order  and  identical  to  the  reference 
trajectory.  Now  if  the  plant  itself  is  not  stable,  i.e. ,  if  the  poly¬ 
nomial  H(z)  has  some  of  its  roots  outside  the  unit  circle,  then  U(z) 
corresponds  to  an  increasing  sequence  u (t) .  This  latter  case  characterizes 
nonminimum  ohase  systems  where  the  cancellation  of  z"*H(z)  in  the  original 
expression  of  ^-|y|  ,  leading  to  (3.10),  is  not  valid  because  H(z)  contains 
unstable  zeros. 

Let  us  come  back  to  the  case  of  imperfect  identification,  H(z)^H(z). 
From  the  transfer  functions  (3.19)-  (3.22)  it  becomes  clear  that  for  the 
system  to  be  stable,  i.e.,  the  output  y(t)  to  be  convergent  and  the  opti¬ 
mum  input  sequence  to  be  bounded,  it  is  necessary  and  sufficient  that: 


Y(z) 
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open  loop:  H(z)  -  az~1H(z)  has  all  its  roots  within  the  unit  circle 

closed  loop:  z"*(l-a)H(z)  +  (l-z~^)H(z)  has  all  its  roots  within 

the  unit  circle 

It  is  easy  to  verify  that  the  above  polynomials  are  characteristic  poly¬ 
nomials  of  the  autoregressive  models  (3.17)  and  (3.18). 

Note  that  in  general  the  exact  plant  transfer  function  H(z)  is  not 

known  and  therefore  one  cannot  evaluate  the  exact  exoression  of  the 
above  polynomial.  However,  if  the  identification  of  H(z)  is  fairly  good, 
then  it  is  expected,  by  the  continuity  theorem  (which  states  that  the 
roots  of  a  polynomial  are  continuous  functions  of  its  coefficients),  that 
the  roots  of: 

H(z)  -  az_1H(z)  (open  loop) 

and  z_1(l-a)H(z)  +  (l-z_1)H(z)  (closed  loop) 

are  close  to  the  roots  of: 

(l-az_1)R(z) 

Therefore,  the  stability  of  the  identified  model  H(z)  implies  the  stability 
of  the  system.  But  when  the  identification  error  becomes  large  so  that 
the  discrepancy  between  H(z)  and  H(z)  becomes  significant,  then  to  deter¬ 
mine  the  stability  of  the  system  one  should  have  recourse  to  robustness 
analysis,  which  is  discussed  in  Section  3.5  and  involves  determining  the  set 
of  plant  polynomials  H(z)  for  which  the  above  characteristic  polynomials 
have  stable  roots. 
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Assume  that  the  uncertainties  of  the  system  are  modeled  as  an  addi¬ 
tive  zero  mean  noise  on  the  plant  output.  That  i s : ^ 

y(t+l)  =  h^u(t)  +  w(t)  (3.25) 

The  prediction  model  and  the  reference  trajectory  being  as  previously 
defined : 

y 

yp(t+l)  =  h^u(t)  open  loop  (3.26) 

yp(t+l)  =  y(t)  +  hT(u(t)  -  u(t-l))  closed  loop  (3.27) 

yp(t+k)  =  aky(t)  +  (l-ak)c  k  =  1 , —  ,T  (3.28) 

The  cost  function  to  minimize  over  a  horizon  of  T  is  now,  instead  of  dp 

E[Jt]  =  e{  ^  (yp(t+j)  -  yr(t+j))2| 

where  E  denotes  the  expectation  operator. 

Following  the  same  line  of  argument  as  in  Section  3.1,  it  is  seen  that 
when  there  is  no  constraint  on  inputs,  the  optimum  input  u*(t)  is  determined 
by  the  equation 

yp(t+i)  =  yr(t+i)  (3.29) 

which  translates  into: 

^Here,  the  output  y(t+l)  may  be  considered  the  measurement  available  to  the 
controller.  It  includes  the  true  plant  output  plus  a  measurement  noise  term. 
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hTu*(t)  =  ay(t)  +  (l-a)c  (open  loop)  (3.30) 

y(t)  +  RT(u*(t)  -  u*(t-l))  =  ay(t)  +  (l-a)c  (closed  loop)  (3.31) 

Similar  to  the  deterministic  case  the  optimum  control  sequence  u*(t) 
is  determined  sequentially  from  the  above  equations,  in  which  the  plant 
output  y(t)  is  observed.  The  difference  with  the  deterministic  case  is 
that  the  present  sequence  of  optimum  control  u*(t)  is  nondeterministic, 
which  results  from  the  noisy  nature  of  the  Diant  output  y(t).  The 
reader  can  easily  verify  that  since  the  noise  has  zero  mean  the  mean 
value  y(t)  of  the  plant  output  and  the  meanu*(t)of  the  optimum  input 
sequence  have  a  deterministic  dynamic  identical  to  the  one  governing  the 
deterministic  case  as  in  Section  3:3.  Thus,  all  the  results  of  the  pre¬ 
vious  section  are  valid  for  the  calculation  of  the  means  u(t)  and  y(t).  It 
remains  to  study  the  variance  behavior,  and  in  particular  the  variance 
of  the  controlled  output  y(t).  Here  it  becomes  necessary  to  differentiate 
between  the  cases  of  white  and  colored  noise. 

3.4.1  White  Noise 

Let  the  noise  be  white,  Gaussian,  and  independent  from  the  input  u(t): 

E[w( t) ]  =  0 

(a2  for  t  =  t 

E[w(t)w(x)]  =  -! 

(O  for  t i  t 

E[w( t) u( t) ]  =  0 

Let  us  first  assume  that  the  identification  of  the  impulse  response  is 
sufficiently  good  to  allow  the  approximation  h^  =  h"*".  Then  the  control 
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equation  becomes: 


yp(t+l)  =  y(t+l)-w(t)  =  ay(t)  +  (l-a)c  (open  loop)  (3.32) 

yp(t+l)  =  y(t)  +  hT(u(t)  -  u(t-l))  =  y(t-l)  - w(t) +w(t-l) 

=  ay(t)  +  (l-a)c  (closed  loop) 

(3.33) 

2  2 

Denoting  by  OqU)  and  ac( t)  the  variances  of  the  open  and  closed  loop 
plant  output  y(t),  one  deduces: 

o2(t+l)  =  a2Og(t)  +  a2  open  loop 
a2(t+l)  =  a2a2(t)  +  2o2(l-a)  closed  loop 

which  result  in  steady  state  variances  of 

open  loop  (3.36) 

closed  loop  (3.37) 

It  is  seen  that  for  the  open  loop  prediction,  the  variance  of  the 

output  decreases  with  a;  that  is,  a  fast  reference  model  (small  a)  results 

in  a  smaller  variance.  But  for  the  closed  loop  case,  the  variance  of  the 
* 

output  is  a  decreasing  function  of  a;  thus,  a  fast  reference  model  results 
in  a  larger  variance.  Figure  9  displays  that  the  breakeven  point  is  for 
a=  V2.  Now, the  practical  implication  of  the  above  results  is  that  when 
the  plant  output  y(t)  is  far  from  the  desired  value  c,  one  can  operate  in 


°o(“) 


°cH 


1-a 

2a! 

1+a 


(3.34) 

(3.35) 
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two  consecutive  phases;  first  use  an  open  loop  prediction  with  fast  ref¬ 
erence  trajectory  (a  small)  and  then  switch  to  a  closed  loop  prediction 
control  with  slow  reference  trajectory.  The  latter  is  necessary,  since 
as  mentioned  before,  the  open  looo  prediction  yields  a  steady  state  bias. 


Now  assume  that  the  identified  impulse  response  hT  differs  from  the 
actual  one  h  .  Then  the  evaluation  of  the  output  variance  involves,  as 
an  intermediate  stage,  the  computation  of  the  covariance  matrix: 

U  (°°)  =  E[u(«)uT(«>)  j 

Or,  equivalently,  the  evaluation  of  E[u(t)u(t-i)]  for  i  =  0,...,N.  It 
can  be  shown  that  this  involves  the  solution  of  a  system  of  linear 

equations  with  unknowns,  this  latter  being  elements  of  U(«>) 

(Anderson  and  Moore,  1979). 

Let  us  consider  the  simple  case  of  constant  gain  mismatch  between  the 
plant  Impulse  response  and  the  model  one,  i.e.. 
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q:  qam 


h 


T 


Then  the  basic  control  equalities  for  open  and  closed  loop  prediction 
become: 


y(t+l)  =  qay(t)  +  q(l-a)c  +  w(t)  open  loop  (3.38) 

y(t+l)  =  [1  -  (l-a)q]y(t)  +  q(l-a)c  +  w(t)  -  w(t-l)  (3.39) 

closed  loop 

The  corresponding  variance  dynamics  are: 

OqU+I)  =  (qa)^OQ(t)  +  a2  open  loop 
a2(t+l)  =  (1  -  (l-a)q)2a2(t)  +  2a2(l-ct)q  closed  loop 

2  2 

whence  the  stability  conditions  for  the  convergence  of  and  oc  are: 

! qot |  <  1  open  loop 
) 2  -  ( l-a)q |  <  l  closed  loop 


(3.40) 

(3.41) 


and  the  corresponding  steady  state  variances  are: 


On  = 


°c  2  -  ( l-a)q 


1  -  (qa)‘ 
2 


2aL 


open  loop 


closed  loop 


From  (3.42),  it  is  clear  that  the  stability  of  both  open  loop  and  closed 
loop  systems  are  guaranteed  for: 
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0  <  q  <  min[ 


1-cx 


that  is, 


9 

q  <  —  for  -1  <  a  <  !/3 

i  t 

q  <  —  for  a  s  '/3 

a 


2  2 

The  plots  of  Figure  10  display  the  variances  oQ  and  as  functions  of  the 
mismatch  gain  q  in  the  above  cases. 


(b)  y3  <  a  <  1 
Figure  10.  Noise  variances  with  mismatch. 


Let  us  first  note  that  the  parameter  a  is  usually  positive,  since  it 
is  very  unlikely  and  undesirable  to  select  oscillatory  reference  trajec¬ 
tories.  Plot  (a)  indicates  that  for  a  fast  reference  model  (0<a<  V3), 
the  open  loop  variance  is  always  smaller  than  the  closed  loop  one  irrespec¬ 
tive  of  the  value  of  the  qain  q  (provided  that  the  system  remains  stable, 


37 


i.e.,  q<2/l-a).  In  contrast,  plot  (b)  shows  that  for  a  slow  reference 
model  (Y3<a<  1),  the  closed  loop  prediction  always  brings  a  smaller 
variance  as  compared  to  the  open  loop  variance,  and  moreover  for  a  gain 
factor  of  less  than  1/a,  the  closed  loop  variance  is  less  than  2o  ,  i.e. 
twice  the  variance  of  the  noise.  Since  in  most  practical  applications 
one  has  to  use  a  closed  loop  scheme  to  assure  the  unbiasedness,  it  is 
desirable  for  the  reference  trajectory  to  have  a  rate  of  growth  a  within 
the  range  of  [Y^,  1).  This  is  not  restrictive  since  in  most  oases,  the 
reference  trajectories  rarely  have  rates  of  growth  less  than  0.6. 

3.4.2  Colored  Noise 

If  the  noise  is  colored,  then  the  relatively  simple  results  of  the 
previous  subsection  are  no  longer  valid.  The  colored  nature  of  the  addi 
tive  noise  on  the  outout  is  often  due  to  the  fact  that  the  plant  output 
y(t)  is  observed  through  a  filter,  usually  of  low  order.  Moreover,  any 
additive  white  process  noise  or  input  disturbance  will  reflect  itself  as 
output  colored  noise. 

To  see  the  type  of  behavior  implied  by  the  colored  character  of  the 
noise,  consider  the  simple  case  of  a  perfectly  identified  system 
(h^=h"^)  in  the  presence  of  a  first  order  Markov  process  noise: 

w(t+l)  =  pw(t)  +  e(t)  |p|  <  1 

p 

e(t)  :  zero  mean  white  noise  of  variance  v 
Tne  stationary  variance  of  the  noise  is  then: 
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The  systems  equations  are  similar  to  the  white  noise  case,  but  the  recursive 
relations  governing  the  variances  are  more  complex. 

For  the  open  loop  prediction  scheme,  it  is  seen  without  much  diffi¬ 
culty  that  the  variance  recursion  is: 

Og(t+l)  =  a2o2(t)  +  o2(t)  +  2a  cov(t) 

where: 

cov ( t ,t)  -  E[(y(t)  -  y(t) )  w( t) ] 
and 

cov(t+l,t+l)  =  aocov(t.t)  +  pa£(t) 

Letting  t-*-°°,  one  deduces  the  stationary  variance  of 

2  _  1+gp  1  2  _  1+ap  1  v2 

0  1-ap  .  2  w  1-ap  .  2  ..  2 

1-g  1-a  1-p 

As  the  plot  of  Figure  11  shows,  for  a  positively  correlated  noise  (p>0), 

2 

the  qualitative  behavior  of  the  variance  Og  with  respect  to  a  remains 
the  same  as  in  the  white  noise  case;  that  is,  an  increase  of  a  (slowing 
the  reference  model)  results  in  a  larger  variance.  But  if  the  noise  is 
negatively  correlated  (p<0),  then  there  is  an  optimum  value  of  a,  i.e., 
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an  optimum  reference  trajectory,  associated  with  each  p.  And  moreover, 

there  is  a  range  of  values  for  a,  resulting  in  an  output  variance  smaller 

2 

than  the  noise  variance  o  . 


Figure  11.  Open  loop  variance,  colored  noise. 


For  the  closed  loop  prediction,  from  the  equation  (3. 33) ,  one  deduces 
the  following  recursions  for  the  variance: 


o^(t+l)  =  a2o^(t)  +  ( l-2p)c^(t)  +  a^(t-l)  +  2acov(t,t)  -  2acov(t,t-l) 


where 


cov  (t,t)  -  E[(y(t)  -  y(t))w(t)] 
cov ( t , t- 1 )  £  F[(y(t)  -  y( t ) )w( t-1 ) ] 


cov(t+l,t+l)  =  apcov(t.t)  +  po2(t)  -  p2a2(t-l) 
cov(t+l,t)  =  acov(t.t)  +  o2(t)  -  pa2(t-l) 

Letting  t  -*■  the  expression  of  the  stationary  variance  is  obtained  as: 

02  =  J_  1-P  2  =  _2__  v2 

1+a  l-p«  w  1+a  1-pa  .  2 

i -i- 

As  the  plots  of  Figure  12  show,  for  a  negative  correlation  factor  p, 

one  has  to  slow  down  the  reference  trajectory  (increase  a)  to  decrease  the 

variance.  For  positive  r,  three  situations  arise  depending  on  the  value 

of  o.  For  o<  '/■>,  the  variance  is  still  a  decreasing  function  of  a.  For 

V3  <  p  <  V2,  the  variance  reaches  its  minimum  for  a  =  ,  and  moreover, 

it  is  always  bounded  by  its  value  at  a=0,  i.e.,  2v  /  1+p.  For  V2<p<l, 

the  situation  is  reversed  in  the  sense  that  the  variance  is  bounded  by  its 

2  2 

value  at  a =  1,  that  is,  by  v  /1-n  .  In  these  two  latter  cases,  slowing 
down  the  reference  model  does  not  necessarily  lead  to  a  smaller  output 
variance  since  the  minimum  is  achieved  for  some  a<l. 


Figure  12.  Closed  loop  variances,  colored  noise. 
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o  2 

Fiqure  13  shows  the  ratio  cT/oq  for  different  values  of  a  and  p.  It  is 
seen  that  the  breakeven  value  of  a,  that  is,  the  value  for  which  the  closed 
loop  and  the  open  loop  variances  are  equal,  moves  toward  the  left  as  the 
noise  correlation  factor  increases.  Therefore,  as  the  noise  becomes 
more  and  more  positively  correlated,  the  range  of  value  of  a  corresponding 
to  a  smaller  closed  loop  variance  becomes  larger. 


The  above  colored  noise  analysis  shows  that  it  is  often  advantageous 
to  Derform  an  output  filtering,  particularly  when  the  output  noise  level 
is  large.  A  proper  choice  of  an  output  filter  decreases  the  output  vari¬ 
ance  while  allowing  the  selection  of  a  faster  reference  trajectory. 

3.5  Robustness 

The  purpose  of  this  section  is  to  formulate  and  study  the  robust¬ 
ness  oroblems  associated  with  the  control  procedures  described  in  the 


previous  sections.  The  control  design  is  said  to  be  robust  if  the  plant 
output  y(t)  converges  to  its  ultimate  desired  value  c,  under  a  range  of 
plant  and  behavior  changes  caused,  for  instance,  by  failure  of  some  of  its 
components  or  by  large  parameter  deviations.  The  reader  should  also 
see  Appendix  B  for  an  application  of  the  ideas  of  this  section  to  a  simple 
missile  example. 

With  regard  to  the  closed  loop  prediction  control  mentioned  previously, 
the  robustness  problem  is  posed  in  terms  of  the  stability  of  the  input 
sequence  u(t)  and  the  plant  output  y(t).  This  is  a  consequence,  as 
mentioned  in  Section  3.3,  of  the  fact  that  the  stability  of  the  closed  loop 
system  guarantees  the  unbiased  convergence  of  the  output  to  the  desired  set 
point  c.  Therefore,  as  long  as  the  structural  changes  of  the  plant,  which 
amount  to  altering  the  impulse  response  [J,  are  such  that 


(z-l)H(z)  +  (l~a)H(z)  (3.43) 

has  all  its  roots  within  the  unit  circle,  the  control  is  robust. 

Let  us  represent  the  relationship  between  the  assumed  impulse  response 
hT  and  the  plant  response  hT  by  a  diagonal  matrix  R: 

hT  =  hTR  or  h.  =  r.h.  i  =  0,...,N-l  (3.44) 


Then  the  characteristic  equation  is  written  as: 


N-l  ~  -i  N-l  _  . 

(z-1)  l  h.z  +  (1-a)  l  r.h.z"1  =  0  (3.45) 

i=0  i=0  11 


Given  the  model  impulse  response  hT,  it  is  of  great  interest  to  determine 


the  range  of  the  matrix  R  such  that  (3.45)  is  stable.  To  do  so  one  has 

N 

to  determine  a  subspace  o  of  IR  ,  such  that  for 
(r0,rr...,rN _})  e  0 

N 

the  polynomial  (3.45)  remains  stable.  The  subspace  6  of  F  determines 

N 

uniquely,  through  (3. 44) ,  another  subspace  6  of  IR  for  impulse  responses 
h  corresponding  to  stable  closed  loon  characteri Stic  polynomials  (3.45). 

This  latter  subspace  is  called  the  domain  of  robustness: 

w  N-l  _  N-l 

$  =  fveffr  |(z-l)  l  h-z  1  +  ( 1  -ct)  l  v.z  1  is  stable) 
i=0  1  i=0  1 

N 

Obviously  the  determination  of  such  a  subspace  of  F  requires  complex  search 
•i  I  :nr  7  tnms ,  and  often,  from  a  practical  point,  an  unrealistic  number  of 
imputations.  But  fortunately,  in  almost  all  practical  problems,  valuable 
additional  information  both  on  the  physical  properties  of  the  plant  and 
on  the  failure  and  perturbation  history  of  its  different  components  is 
available.  The  appropriate  use  of  this  information  reduces  substc  ly 
the  size  and  dimension  of  the  search  space.  For  instance,  if  the  robust¬ 
ness  aoainst  identification  errors  or  small  perturbations  is  considered, 
then  the  coefficients  r.  are  of  the  form: 

r .  *  1  +  e .  i  =  0,...  ,N-1 

where  ci  are  small  and  their  upper  bounds  can  be  estimated  from  the  particular 
identification  scheme  in  ooeration.  If  the  robustness  against  failure  of 
some  comoonents--sensors  or  activators--is  of  concern,  as  in  Ackermann  (1979), 
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then  only  a  few  elements  of  the  impulse  response  change  drastically;  that 
is,  most  of  the  coefficients  r..  are  equal  to  1.  In  addition  to  this  re¬ 
duction  of  the  search  space,  various  approximation  procedures  can  be  used 
which  result  generally  in  more  conservative  estimates  of  the  robustness 
domain  <p  but  have  the  advantage  of  computational  efficiency. 

Consider  the  case  where  the  mismatch  between  the  plant  impulse  re¬ 
sponse  and  the  model  is  limited  to  pure  gain,  that  is: 

R  =  ql  or  r.  =  q,  i  =  0,. . . ,N-1 

The  characteristic  polynomial  of  the  closed  loop  prediction  scheme  be¬ 
comes  : 


(z-l)H(z)  +  (l-a)qH(z) 

Since  the  model  H ( z )  is  stable,  the  stabi 1 ity--i.e. ,  robustness--condition 
becomes : 


1  -  q(l-ot)  |  <  1 


that  is. 


However,  in  practice  one  would  limit  the  range  of  admissible  gain  q  to 
1  12 

(0,y- -),  since  the  range  {py  ,  py)  results  in  an  oscillatory  response 
for  y(t) ,  which  is  not  desirable.  A  similar  analysis  can  be  carried  out  to 
determine  the  phase  margin  for  the  controller. 
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3.6  Robustness  Performance  Index 

So  far  we  have  been  concerned  with  robustness  in  an  absolute  sense, 
that  is,  from  the  viewpoint  of  asymptotic  convergence  of  the  plant  out¬ 
put  to  its  desired  value.  Such  a  viewpoint  does  not  provide  an  assess¬ 
ment  of  the  relative  robustness  performance  of  different  control  schemes 
under  various  conditions.  In  particular,  the  effects  of  perturbations 
and  components  failure  on  the  speed  of  convergence  are  not  apparent  at 
all,  and  clearly  this  is  very  important  for  industrial  applications. 

In  light  of  the  above  considerations,  it  becomes  necessary  to  define 
some  type  of  "measure"  for  the  performance  of  the  control  scheme  with 
regard  to  robustness  (Mehra  et  aK ,  1979).  Here  we  propose  two  possible 
performance  indices.  First  consider: 

D(y(t),c|hT=  hT,a) 

P'=  - - t  »  (3.46) 

Sup  D(y(t),c|h  R=  h  ,a) 

Ree 

where  D  denotes  the  total  weighted  distance  between  the  plant  output  y(t) 
and  the  desired  set  point  c: 

00  ? 

D  =  l  (y(t)  -  z)  wt  wt  >  0 
t=0  1  1 

A  decreasing  sequence  of  weight  w^  would  enhance  the  importance  of  an  early 
convergence  to  c.  The  numerator  of  the  index  is  the  distance  conditioned 
on  the  perfect  knowledge  of  the  plant  (R  =  h)  and  the  use  of  a  reference 
trajectory  with  rate  a.  In  the  denominator,  first  the  distance  is  computed 
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assuming  that  the  discrepancy  between  the  plant  h  and  the  model  h  can  be 
represented  by  h  =  hR,  and  the  reference  trajectory  is  specified  by  a. 

Then  the  maximum  of  this  distance  caused  by  variation  of  R  within  the 
range  9  is  considered.  Both  numerator  and  denominator  depend  on  the  ref¬ 
erence  trajectory  through  the  value  of  a.  Obviously  the  value  of  the 
performance  index  is  within  the  range  (0,1),  and  the  performance  improves 
with  increasing  p'. 

It  is  possible,  and  desirable,  to  define  a  performance  index  inde¬ 
pendent  from  the  particular  reference  trajectory.  Such  an  index  can  be 
used  in  comparison  of  MAC  with  other  control  schemes,  such  as  LQ.  One 
can  define  such  an  index  by: 

Inf  D(y(t) ,c,a|hT  =  h) 

pi.  _  aeA _ ~ 

Sup  Inf  D(y(t),c,a|h^R=  hT) 

Re 6  aeA 

where  ac  (0,1)  is  the  range  of  admissible  reference  trajectories.  The 
main  difference  with  the  previous  case  is  that  the  controller  is  free  to 
choose  the  best  reference  trajectory  in  the  sense  of  minimizing  the  dis¬ 
tance  between  the  plant  output  and  the  desired  path  c.  Here  the  numerator 
denotes  the  minimum  distance  achieved  by  a  proper  choice  of  a,  under  the 
condition  of  perfect  identifiabil ity.  In  the  denominator,  for  each  choice 
of  R,  the  minimum  distance  is  achieved  by  a  proper  choice  of  a;  then  the 
supremum  is  taken  when  R  varies  in  the  range  9.  ^ 

The  above  performance  Indices  are  defined  for  deterministic  models, 
but  their  generalization  to  the  case  corresponding  to  additive  output 
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noise  is  straightforward  by  replacing  the  distances  D  with  their  expec¬ 
tations.  As  an  example,  consider  the  case  of  pure  gain  mismatch  q  and 
additive  white  noise,  then  the  indices  can  be  easily  computed  to  have 
the  following  expressions: 


=  1nf 


e  =  c  £°’I^ 


Sup  Inf 
qe0  ae^ 


|  0  (q1,q2 ^ c  i — ct ^ 

(  A  =  (a14a?) c  (0,1) 


To  see  the  behavior  of  p  with  respect  to  a,  consider: 


dp'  _  _d_  P 

da  da 


2-(l-a)q2]  _  2(q2-l) 
"  (1+a)2 


For  q2  >  1 ,  that  is,  if  the  maximum  possible  mismatch  gain  is  larger  than 
unity  (which  is  usually  the  case),  then  slowing  down  the  reference  trajec¬ 
tory  (decreasing  a)  improves  the  robustness  of  the  system. 

The  expression  of  p"  involves  the  maximum  value  of  the  gain  and  also 
the  maximum  value  of  a  which  corresponds  to  the  slowest  admissible  ref¬ 
erence  trajectory: 

_  2-(l-a2)q2 


which  decreases  linearly  when  q,  varies  from  1  to  ‘  ■  . 

It  is  interesting  to  note  that  the  above  results  on  robustness  en- 
hancemer»t  via  slowing  down  the  reference  trajectory  have  been  confirmed 
in  practical  applications. 
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In  most  applications,  the  input  u(t)  is  not  free  of  constraints. 

In  general  there  are  both  amplitude  and  rate  constraints  imposed  by 
technical  and  cost  considerations.  Under  such  constraints,  the  results 
of  Section  3.2  must  be  revised;  in  particular,  the  minimization  of  the 
distance  Jy  over  a  horizon  of  length  T  (T  >  1)  does  not  necessarily  yield 
the  same  optimal  input  u*(t)  as  the  minimization  of  Jj  (T =  1) . 

It  can  be  shown  that,  under  rather  weak  conditions,  the  stability 
of  the  system  free  of  input  constraints  implies  the  stability  of  the 
input  constrained  system.  In  fact,  Praly  (1979)  (see  Appendix  A)  has 
proved  a  more  general  result  pertaining  to  nonlinearity  in  the  recursive 
control  equation  determining  u(t).  Let  us  describe  briefly  his  result, 
a  complete  proof  of  which  involves  mathematical  technicalities  and  can  be 
found  in  Appendix  A.  When  there  are  no  constraints,  the  recursive  equa¬ 
tion  generating  the  inputs  u(t)  is  linear,  of  the  type: 

i  r  N-i  *| 

U*(t+1)  =  r—  !-  I  h,u*(t-j)  +  ay(t-l)  +  (l-a)c  (3.47) 

0  L  j=l  J  J 

Assume  that  the  above  recursion  is  stable  and  denote  its  steady  state 
solution  by  u*.  Now  consider  a  nonlinear  time  varying  function  f  (•) 
and  assume  that  instead  of  (3.47)  the  recursion  governing  the  input 


sequence  is: 


Then  in  order  for  the  stability  of  (3.47)  to  imply  the  stability  of 
(3.48)  it- is  necessary  that  at  every  time  t,  and  for  any  v: 


|ft(u*  +  v)  -  U*|  <  V [ 


(3.49) 


where  r  is  the  greatest  modulus  of  the  roots  of  the  characteristic  equation 
corresponding  to  (3.47). 

In  the  case  of  constant  amnlitude  bounds  on  the  input  u*(t),  the 
functions  ft(*)  are  time  independent  and  involve  saturation  type  of  non- 
linearities: 


f(v)  ^ 
with  u* 


M  for  v  >  M 
•  v  for  m  <  v  <  M 
m  for  v  <  m 

e  [m,M] 


It  can  be  shown  that  such  functions  verify  the  inequality  (3.49)  when  the 
recursion  (3.47) is  stable.  Therefore,  the  nonlinearities  introduced  by 
amplitude  constraints  do  not  make  the  system  unstable. 

The  above  result  is  important  since  it  reduces  the  absolute  stability 
analysis  of  a  MAC  driven  system  with  amplitude  constraints  to  the  one  of  the 
unconstrained  system,  which  is  relatively  simple,  thanks  to  its  linearity. 
It  should  be  noted  that  the  necessary  condition  (3.49)  is  relatively  gen¬ 
eral.  A  rather  large  class  of  nonlinear  time  varying  functions  verify 
such  an  inequality.  These  functions  may  reflect  different  types  of 
constraints  (such  as  time  varying  rate  and  amplitude  constraints),  and 
therefore  the  inequality  (3.49)  can  be  used  to  relate  the  stability  analysis 
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of  a  system  with  more  general  types  of  constraints  to  that  of  the  un¬ 
constrained  system. 


3.8  Summary 

A  mathematical  framework  for  the  analysis  of  an  idealized  MAC 
scheme  has  been  presented.  The  main  operations  performed  by  MAC  and  the 
corresponding  components  have  been  described.  For  single  input-single 
output  systems,  specific  results  pertaining  to  the  stability  and  robust¬ 
ness  have  been  derived.  An  example  of  stability  analysis  of  a  MAC 
control  system  for  the  single  input-single  output  case  is  presented  in 
Appendix  A. 


SECTION  IV 


MAC  FOR  NONMINIMUM  PHASE  AND  TIME  DELAY  SYSTEMS 


Overview 

In  this  section,  the  control  of  time  delay  and  nonminimum  phase 
systems  is  discussed  within  the  framework  of  MAC.  It  is  shown  that  for 
the  control  of  such  systems,  MAC  offers  a  relatively  wide  range  of  solu¬ 
tions  with  various  levels  of  performance.  The  implementations  of  these 
solutions  are  fairly  simple  and  often  involve  minor  alteration  of  the 
existing  MAC  software.  However,  from  a  theoretical  point  of  view,  there 
is  a  major  difference  between  minimum  phase  and  nonminimum  phase  systems, 
namely,  the  inverse  of  the  former  is  stable  whereas  the  inverse  of  the 
latter  is  unstable.  This  necessitates  the  use  of  different  models  for 
control  computations  and  for  prediction  in  MAC.  The  control  model  has  the 
property  of  possessing  a  stable  inverse  whereas  the  prediction  model  must 
be  as  close  as  possible  to  the  true  model.  In  this  section,  an  optimal 
control  model  for  use  in  MAC  will  be  derived  and  its  performance  with  other 
commonly  used  techniques  for  nonminimum  phase  systems  will  be  compared. 

The  organization  of  the  present  section  is  as  follows.  Section  4.1 
presents  a  historical  background  and  4.2  presents  a  brief  review  of  MAC, 
describing  its  essential  features  for  the  single  input-single  output  case. 
Section  4.3  discusses  the  general  nature  of  nonminimum  phase  systems  in 
the  framework  of  MAC.  It  is  also  concerned  with  qualitative  and  conceptual 
descriptions  of  the  approaches  of  Sections  4.4  and  4.5.  Section  4.4 


describes  the  "least  squares  solution"  for  nonminimum  phase  systems.  For 
the  sake  of  clarity  the  detailed  mathematical  derivation  of  this  latter 
solution  is  given  in  Appendix  C.  Section  4.5  presents  a  minimum  phase 
solution  based  on  removal  of  unstable  poles.  This  is  basically  a  pole 
placement  type  of  solution.  Section  4.6  describes  more  heuristic  approaches 
which  are  nonetheless  of  great  use  and  importance  for  applications.  These 
heuristic  approaches  are  basically  trial  and  error  types  of  solutions 
involving  both  cost  coefficients  and  prediction  interval  adjustments. 

4 . 1  Background 

Control  theory,  from  its  very  outset,  has  paid  particular  attention 
to  the  control  of  nonminimum  phase  systems.  This  interest  has  been  moti¬ 
vated  by  the  relatively  frequent  occurrence  of  nonminimum  phase  properties 
in  various  engineering  control  problems,  particularly  with  multivariate 
systems . 

All  the  available  optimum  control  schemes--! .e. ,  Linear  Quadratic 
control  (LQ),  Minimum  Variance  control  (MV),  Self-Tuning  (ST)  and  Model 
Adaptive  control  (MA)--have  developed  their  own  methodologies  to  tackle 
nonminimum  phase  systems  (Astrom,  1970;  Peterka,  1972;  Astrom  and  Witten- 
mark,  1974;  Ogata,  1970;  Kwaakernak  and  Sivan,  1972).  Despite  the  fact 
that  these  methodologies  stem  from  different  control  philosophies,  they 
have  some  common  features.  In  particular,  all  of  them  lead  in  one  way  or 
another  to  suboptimal  solutions  with  bounded  inputs.  The  optimality  of 
the  solution  is  usually  traded  off  for  a  meaningful  realizable  suboptimal 
solution  (Astrom,  1970;  Peterka,  1972). 
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The  present  section  treats  nonminimum  phase  systems  in  the  framework 
of  MAC.  As  will  be  seen,  the  MAC  methodologies  for  such  systems  offer  a 


set  of  possible  solutions  ranging  from  the  "best  suboptimal  control," 
which  involves  an  off-line  solution  of  a  Riccati  equation,  to  a  more 
heuristic  type  of  control  involving  pole  placement  and/or  cost  coefficient 
adjustment.  The  general  properties  of  these  solutions  are  similar  in 
many  respects  to  the  ones  derived  from  previously  mentioned  control  schemes 
(LQ,  MV,  ST,  MA).  In  particular  it  is  shown  that  there  is  always  a  trade¬ 
off  between  the  optimality  of  a  solution  and  the  boundedness  of  its  cor¬ 
responding  input.  The  practical  implementation  of  the  MAC  solutions  for 
nonminimum  phase  systems  often  needs  only  minor  alterations  of  the  existing 
MAC  software  and  hence  may  be  realized  at  relatively  low  cost. 

4.2  Brief  Review  of  MAC 

We  recall  from  Section  III  that  in  MAC  the  true  plant  and  its  mathe¬ 
matical  model  are  represented  by  their  impulse  responses  hT  and  hT: 

y(t+l)  =  hTu(t)  Plant 

(4.1) 

yM(t+l)  =  hTu{t)  Model 

N  ^  N  N 

with  heF  ,  heF  ,  u(t)  eF  .  The  vector  u(t)  represents  the  past  N  inputs 
u(t), . . . ,u(t-N+l) . 

Based  on  this  mathematical  model  and  on  the  observed  output  of  the 
plant,  a  linear  prediction  scheme  is  conceived: 
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(4.2) 


yp(t+l)  =  yM(t+l)  +  [y(t)-yM(t)] 

=  y(t)  +  hT[u(t)  -  u(t-l)] 

The  controller  also  generates  a  reference  trajectory  y^,  which  is 
initialized  on  the  observed  output  of  the  plant: 

yr(t+k)  =  aky(t)  +  (l-ak)c;  |a|  <1 

where  c  is  the  ultimate  desired  value  to  which  the  plant  output  must 
converge. 

The  control  strategy  of  MAC  is  to  force  the  predicted  output  yp 
to  follow,  as  closely  as  possible,  the  reference  trajectory  yr(t+k)  over 
T  future  periods  (k  =  l,...,T).  Hence  the  controller  must  determine  a  set 
of  T  future  inputs  u(t), . . . ,u(t+T-l)  such  that  the  weighted  distance: 

T  2 

JT  =  ,nyp(t+i)-yr(t+i)]  w.;  w^  >  0  (4.3) 

is  minimized.  It  was  shown  in  Section  III  that  for  a  given  h,  R,  a  and 
c  there  exists  a  unique  sequence  of  inputs  such  that  the  above  minimum  is 
zero,  i.e.,  the  predicted  output  yp  follows  exactly  the  reference  trajec¬ 
tory.  In  particular,  if  fi  =  h,  the  corresponding  optimum  input  sequence 
forces  the  plant  output  y(t)  to  follow  exactly  the  reference  trajectory 
yr.  Therefore,  if  h  is  known,  the  best  choice  for  h,  in  the  sense  of 
minimizing  the  distance  between  the  plant  output  and  the  reference  tra¬ 
jectory,  is  h. 

For  a  one  step  minimization  (T  =  1),  the  optimum  input  u*(t),  mini¬ 
mizing  Jj,  is  determined  by: 
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4 . 3  Nonminimum-Phase  System--Statement  of  the  Problem 

Let  us  assume  that  the  impulse  response  H(z)  of  the  plant  is 
perfectly  known.  Then,  equating  the  model's  response  H(z)  with  the 
plant's  response,  we  get  from  equations  (4.5)  and  (4.6), 


Y(z)  = 

H(z) 

(z-u)  HTzT 

(l-O 

C  ( z )  =  £-C(z) 

(4.7) 

U(z)  = 

rz-aTHTzy 

(1-0 

C(z) 

(4.8) 

Y(z) 

Note  that  the  Diant  output  tracking  transfer  (i.e.,  )  is  of  first 

order  with  unit  gain  and  a  single  pole  at  i.  This  means  that  the  plant 
outout  follows  exactly  the  reference  trajectory.  The  plant  impulse 
response  H(z)  does  not  intervene  in  the  expression  of  Y(z).  But  its 
inverse  intervenes  in  the  expression  of  the  plant  input  U(z).  The  zeros 
cf  H(z)  become  poles  of  the  plant  input  tracking  transfer  (i.e.,  )• 

If  H ( z )  has  some  of  its  zeros  outside  the  unit  circle,  that  is,  if 
the  system  is  nonminimum  phase,  it  is  clear  from  equation  (4.8)  that  the 
plant  input  u(t)  becomes  unbounded.  However,  the  output  y(t)  remains 
bounded  and  follows  the  reference  path  as  long  as  none  of  its  unstable 
poles  (corresponding  to  the  cancelled  H(z)  of  equation  (4.?))  are  ex¬ 
cited.  It  should  be  noted  that  practical  considerations  such  as  physi¬ 
cal  constraints,  cost,  etc.  will  always  impose  an  upper  bound  on  the 
norm  of  u(t).  Once  this  bound  is  reached,  the  cancellation  of  H(z)  in 
equation  (4.7)  is  no  longer  perfect  and  the  output  y(t)  starts  to  di¬ 
verge  from  the  reference  path.  Figure  15  displays  such  a  situation. 
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When  H(z)  is  not  known  precisely,  which  is  usually  the  case,  the 
situation  is  worse  in  the  sense  that  both  y(t)  and  u(t)  diverge.  This 
is  seen  from  the  expression  for  Y (2 )  and  U(z)  in  equations  (4.5)  and 
(4.6).  The  denominator  [(l-a)H(z)  +  (z-l)fl(z)]  has  its  zeros  close  to  the 
zeros  of  (z-a)H(z)  (since  H(z)  is  close  to  H(z)).  Therefore,  Y ( z )  is 
unstable  and  divergent.  (For  further  discussion  of  nonminimum  phase 
systems,  see  Astrom,  1970;  Peterka,  1972;  Ogata,  1970;  Astrom  and  Witten- 
mark,  1974.) 

From  the  above  analysis  we  deduce  that  unlike  the  minimum  phase 
case  where  the  natural  choice  of  h^  (or  equivalently  H(z))  is  a  vector 
identical  or  close  to  hT,  for  the  nonminimum  phase  case  this  choice  of 
hT  does  not  lead  to  a  stable  and  realizable  control  strategy.  In  fact, 
in  order  to  have  a  stable  control  strategy  which  leads  y(t)  to  its  desired 
final  value  c,  it  is  necessary  to  choose  an  H(z)  such  that  the  polynomial 
[H(z)(l-a)  +  (z-a)H(z)]  has  all  its  roots  within  the  unit  circle.  This 
results  in  a  stable  input  sequence  u(t),  which  tends  toward  a  steady 
state  causing  the  term  h^[u(t)  -u(t-l)]  of  equation  (4.4)  to  converge 
to  zero.  Then  as  t -*•»,  lim  yp(t+l)  =  lim  y(t)  =  lim  yf(t+l)  =  c.  Therefore, 
a  proper  choice  of  H(z)  (i.e.,  h"*")  will  ensure  the  ultimate  convergence 
of  the  output  y(t)  to  the  desired  value  c  while  keeping  the  input  sequence 
bounded.  However,  as  one  might  expect,  the  ultimate  convergence  of 
y(t)  to  the  desired  value  c  is  not  done  at  free  cost;  since  ^  h  ,  the 
transitional  behavior  of  y(t)  is  different  from  that  of  Yp(t)  (and  hence 
from  the  reference  model  yf(t)  which  always  equals  yp(t)  under  the  present 
control  strategy).  It  is  worth  noting  that  by  choosing  h^  equal  to  h^. 
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the  inverse  situation  occurs;  that  is,  the  output  y(t)  follows  perfectly 
yp(t)  (and  hence  the  reference  model  yf)  during  the  transitional  phase, 
but  diverges  from  it  as  t  increases  (Figure  15). 

Obviously,  there  are  an  infinite  number  of  vectors  hT  resulting 
in  a  stable  characteristic  polynomial  [H(z) (1-a)  +  (z-l)H(z)] .  Among 
all  the  admissible  hT  one  may  choose  those  which  satisfy  some  secondary 
criteria.  A  natural  criterion  is,  for  example,  the  minimization  of  the 
total  distance  between  the  actual  output  y(t)  and  the  reference  model 
yr(t): 


min  lCy(j)  -yr(j)]2 
j=l  r 

Other  criteria  may  involve  direct  pole  placement  for  the  characteristic 
polynomial  [H(z) (1-a)  + H(z) (z-1)],  or  equivalently  the  design  of  the 
dynamic  of  the  error: 

e(t)  =  yr(t)  -y(t) 

Sections  4.4  and  4.5  discuss  the  above  cases. 

Before  ending  this  section  v/e  note  that  the  use  of  an  hT  in  equation 
(4.7)  which  is  consciously  different  from  h^  implies  that  the  mathematical 
model  yM  is  no  longer  attempting  to  "represent"  the  system.  Therefore, 
instead  of  y^  we  introduce  a  new  variable  e(t)  which  measures  the  error 
between  the  plant  output  y(t+l)  and  the  reference  output  yr(t+l): 
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e(t+l)  =  yr(t+l)  -  y(t+i) 

=  yp(t+l)  -  y(t+i) 

=  [hT-hT][u(t)-u(t-l)] 

=  d[u(t)  -  u(t-l)] 

with  d  =  hT-h*\  As  we  can  see,  the  vector  h"1"  directly  influences  the 
dynamic  of  the  error. 

4.4  Least  Square  Solution 

The  present  section  deals  with  the  determination  of  a  vector  h^ 
such  that  the  total  Euclidean  distance  between  the  actual  output  y(t)  and 
the  reference  model  y^  is  minimized.  Minimize: 

J  ■  >s  l  IrU) -  >»  l  c2(J)  (4.9) 

j=i  r  j-i 

-T 

Hence  the  problem  can  be  stated  as  follows:  determine  a  vector  h 
such  that  the  characteristic  equation  [ ( l-a)H(z)  +  (z-l)H(z)]  has  all  its 
roots  within  the  cnit  circle  and  J  is  a  minimum. 

To  solve  the  above  problem,  as  one  may  expect,  we  will  use  the  jargon 
of  linear  quadratic  control.  But  first  we  have  to  formulate  the  problem 
in  a  framework  suitable  for  a  linear  quadratic  approach.  For  the  sake 
of  simplicity,  and  without  loss  of  generality,  we  may  assume  that  c =  0, 
i.e.,  the  desired  final  value  of  the  output  is  zero.  Then  the  equations 
which  determine  the  optimum  sequence  of  input  u(t)  (in  the  sense  of 
matching  yp(t+l)  with  yf( t+1))  are,  by  equation  (4.4): 
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yr(t+l)  =  ay(t) +  hT(u(t) -u(t-l))  =  yp( t+1) 

In  terms  of  the  sequence  u(t),  we  have: 

i^uft)  =  [hj  -  h2  -  (l-a)h1]u(t-l)  +  . . .  +  [h^  -  KN_2  -  (I-o)hN_1]u(t-N) 
+  [hN-  <l-ct)hN]u(t-N-l)  (4.10) 

The  expression  of  the  error  e(t)  in  terms  of  u(t)  is  derived  as 
follows: 

e(t)  =  y(t)-yr(t)  =  y(t)  -  ay(t-l)  =  hTu(t-l)  -  ahTu(t-2) 

=  Chj»  h2  -  ah^,  ...»  ”®^N^  ^(T~l) 

u(t-N-l) 

u(t-N-2)_ 

or  equivalently, 

e(t)  =  gT  "  u(t-l) 
u(t-N-l 

Note  that  the  dimension  of  g^  is  (N+l),  when  the  dimension  of  h^  is  N. 

Let  us  introduce  the  (N+l)-dimensional  vector  §(t)  defined  as: 

§T(t)  =  [uT(t)  {  u(t-N-l)]  =  [u(t),  u(t-l) , . . . ,  u(t-N),  u(t-N-l)] 
Then  e(t)  and  J  can  be  written  as 
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with  g  -  Chj,  h2  -  oihj,  •••»  h^  ■*  cih|^_j»  -(xh|^] 


(4.11) 


i 


(4.12) 


e(t)  =  gTC(t) 

J  =  lim  CJs§T(t)qqT?(t)  +  £  CT( j)qgT§(j)] 
t-*»  j=l 


Equation  (4.10)  can  be  written  in  terms  of  §(t)  as: 

§(t)  =  U§(t-1)  +b1{[LT  j  0]|(t-l)}  (4.13) 

with  U,  bj,  and  LT  defined  as: 

(N+l)  (4.14) 

L  =  ~  [hj  +  ^2“ (l-a)hj,. .. ,  h^_ j  +  h^-( 1-a) ^ »  h^- ( l-a)h^] 
hi 


o' 

o 

o 

o, 

> 

T 

10  0  0 

(N+l) 

0 

0  10  0 

0 

0  0  10 

bl  ■ 

0 

0 

0 


Equation  (4.13)  is  a  state  equation  with  a  classical  state  feedback 
control  law.  It  is  clear  that  the  determination  of  the  optimum  hT  which 
minimizes  J  is  equivalent  to  the  determination  of  a  linear  feedback  gain 
[L^jO]  of  equation  (4.13)  minimizing  J.  At  this  point  there  is  sufficient 
motivation  to  formulate  the  problem  in  the  following  way. 

Find  the  optimum  input  sequence  u*(t)  which  minimizes 


k-1 


°  ‘  1k^'iC5T(k>qqT5(k)+i15T(Jl<l'lT5(J)3 


(4.15a) 
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where  the  state  equation  is: 


§(t)  =  Ug(t-l)  +  b1u*(t-l) 


It  is  easy  to  check  that  the  above  system  is  controllable.  The 
optimum  control  u*(t)  is  a  linear  function  of  the  state  vector  £(t)  :  u 
LT§( t )  (Sage  and  White,  1977;  Murata,  1977).  It  remains  to  show  that 
L  has  the  property  of  having  its  last  element  L^+1  equal  to  zero.  The 
following  lemma  establishes  that  result. 

Lemma:  The  optimal  input  sequence  u*(t)  for  equation  (4.15)  does 
not  depend  on  the  last  element  of  the  state  vector  f;(t). 

Proof:  By  linear  quadratic  control  theory,  we  know  that  u*(t)  has 
the  following  expression  in  terms  of  C(t): 


u*(t)  =  -{[R  +  b1TPb1]'1b1TPU}g(t)  =  L§(t) 


where  R  is  the  cost  matrix  of  the  input  u(t)  (here  R  =  0  since  J  depends 
only  on  the  state  vector  and  P  is  the  steady  state  solution  of  a  matrix 
Riccati  equation).  Then, 


bl  Pbl  =  P11 


b^PU  =  [Pjj,  P12»  •  Pi(N+l)JU  =  ^P12’  P13*  •••’  P 


and  hence 


1(N+1)  ’ 


Let  us  recapitulate  the  results  so  far.  It  has  been  shown  that  the 
determination  of  the  best  vector  h"*"  in  the  sense  of  minimizing  the  total 
Euclidean  distance  between  the  plant  output  and  the  reference  trajectory 
while  keeping  the  control  sequence  u*(t)  bounded  is  equivalent  to  the 

A 

determination  of  the  gain  matrix  L  of  the  linear  quadratic  control  prob¬ 
lem  (4.15).  The  determination  of  L  involves  the  solution  of  a  Riccati 

/s  ~T 

equation  of  order  N+l.  Once  L  is  known  h  can  be  computed  uniquely  from 
the  last  equation  of  (4.14).  It  should  be  noted  that  h^  depends  uniquely 
on  h^  and  «  (the  rate  of  growth  of  the  reference  model).  Since  the  de¬ 
termination  of  hT  is  done  off-line,  even  for  large  values  of  N  (of  order 

~  T 

80)  the  vector  h  can  be  computed  without  excessive  cost.  The  linear 
quadratic  control  problem  of  (4.15)  is  further  detailed  in  Appendix  C. 

We  end  the  present  section  by  applying  the  above  considerations  to  the 
example  of  the  previous  section  (see  Figure  16). 

4 . 5  Direct  Pole  Placement 

In  addition  to  the  above  optimal  choice  in  the  sense  of  quadratic 
error  minimization,  there  are  an  infinity  of  other  possible  choices  of 
h^  leading  to  stable  control  solutions.  The  vector  h^  is  to  be  chosen 
such  that  the  polynomial  [H(z)(l-a)  +  (z-l)H(z)]  has  all  its  roots  confined 
within  the  unit  circle.  Therefore,  determining  H(z)  amounts  to  identi¬ 
fying  the  above  characteristic  polynomial  with  an  arbitrary  stable 
polynomial  of  order  N+l. 

T  ~T 

However,  going  back  to  the  original  definition  of  h  and  h  ,  and 

~T  .  T 

noting  that  for  minimum  phase  systems  one  equates  h  with  h  ,  it  becomes 
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natural  to  conceive  a  design  procedure  for  H(z)  which  remains  consistent 
with  the  minimum  phase  case.  Loosely  speaking,  H(z)  should  retain  "as 
much  as  possible"  of  the  dynamic  properties  of  H(z),  while  stabilizing 
the  characteristic  polynomial  [H(z) (1-a)  + (z-l)H(z)] .  One  way  of  achieving 
that  objective  is  to  maintain  the  stable  roots  of  H(z)  and  to  alter  only 
the  unstable  ones.  Let  us  write  H(z)  as: 

H(z)  =  B0s(z)BQu(z)  (4.16) 

where  B0s(z)  contains  all  the  stable  roots  of  H(z)  and  BQu(z)  the 
remaining  unstable  ones.  Then  let 

fdeg  Bn_  =  N-£ 

H(z)  =  Bfl  (z)Bm.(z)  0S  (4.17) 

ldeg  BMs  -  l 

where  BQs(z)  has  the  same  stable  zeros  as  H(z).  The  characteristic  poly¬ 
nomial  [H(z)(l-a)  +  (z-l)H(z)]  becomes 

B0s(z)[B0u(z)(1"a)  +  (4.18) 

BMs ( z )  must  be  cbosen  so  that  polynomial 

B0u(z)(l-a)  +  (z-l)BMs(z)  (4.19) 

is  stable.  Once  B^s(z)  is  specified,  the  polynomial  H(z)  and  thus  the 
vector  hT  are  determined.  The  resulting  control  sequence  u(t)  will  display 
' ' '  'he  stable  modes  of  H(z)  in  addition  to  the  one  of  [BQu(z)(l-a)  + (z-l)BMs(z)]. 
''owe  for  this  latter  polynomial  is  the  one  having  (£+1)  zeros  at 
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the  origin: 


i 


* 

* 


B0u(z)(l-a)  +  (z-1)BMs(z)  =  CjZU1 
or  (1-z)BMs(z)  =  C1z£+1  -  (l-a)BQs(z) 


(4.20) 


with  being  a  real  constant. 

The  above  choice  of  BMs(z)  will  lead  to  a  fast  convergence  of  y(t), 
although  not  a  very  smooth  one.  Other  choices  of  (4.19)  may  lead  to  a 
smoother  convergence  of  y(t).  In  particular,  BMs(z)  may  be  determined 
so  that  the  polynomial  (4.19)  becomes  identical  to: 

C2zk(z£''1"k  -  r) 

with  being  a  real  constant,  G<k<£+1,  and  |r|  «  1.  That  is,  the  (£+1) 
unstable  poles  of  the  original  system  are  replaced  by  k  poles  at  the  origin 
and  (2.+l-k)  poles  of  module  |r|1^+1_k.  Figures  17  and  18  display  the 
application  of  the  above  considerations  to  the  example  of  Sections  4.4 
and  4.5. 

Before  ending  this  section,  let  us  notice  that:  (i)  the  polynomial 
Bq^z)  may  also  include  stable  poles  which  are  close  to  the  unit  circle. 

By  doing  so,  one  is  assured  that  "small"  perturbations  of  H(z)  will  not 
destabilize  the  controlled  system. 

(ii)  The  main  advantage  of  the  "direct  pole  placement"  approach 
over  the  linear  quadratic  one  lies  in  its  greater  flexibility  as  compared 
with  the  rigidity  of  the  unique  aspect  of  the  LQ  solution.  This  flexi¬ 
bility,  which  is  the  result  of  the  relative  freedom  involved  in  the  choice 
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4>  y(t)-yr(t) 


hT  =  [l.t  -.9,  -1.26,  1.806,  -1.1092,  .2856],  hT  =  [-.096,  -1.904, 
.1915,  -.1685,  .0605],  Poles  of  the  nonminimum  phase  system: 
1.2,  1.4.  Poles  of  [H(z)(l-a)  +  (z-l)H(z)]:  .3±5j,  .5,  0,  0,  0 

a  =  .8,  c  =  0. 


Figure  17.  Stabilization  by  pole  placement. 


of  BMs(z),  allows  one  to  influence  some  important  aspects  of  the  converging 
output  y(t)  such  as  its  smoothness,  overshoot  and  undesired  oscillations. 

In  many  practical  problems  one  is  willing  to  sacrifice  the  "minimum 
quadratic  error"  nature  of  the  linear  quadratic  solution  for  an  improve¬ 
ment  in  the  other  more  important  features,  such  as  those  mentioned  above. 

4.6  Constant  Input  Over  p  Steps 

In  the  present  section  we  describe  another  somewhat  more  heuristic 
approach  to  the  control  of  nonminimum  phase  systems.  The  present  approach 
is  different  from  the  above  in  that  it  requires  some  alteration  of  the 
structure  of  MAC  strategy.  However,  at  the  level  of  implementation  these 
necessary  alterations  are  done  via  minor  changes  in  the  existing  MAC  soft¬ 
ware. 

Let  us  recall  that  in  the  regular  MAC,  at  each  period  t,  one  looks 
for  p  future  optimum  inputs  |u(t),  u(t+l),...,  u(t+p-l)}  such  that  the 
sum 


l  w.|y  (t+i-1)  - oc1y(t ) | 2  =  l  wi|y(t+i) +h(u(t+i) 
i=0  1  r  i=0  1 

-  u ( t+i - 1 ) )  -  aiy(t)|2. 


Wj  >  0 

(4.21) 


is  minimized. 

Now  let  us  impose  the  restriction  of  constant  control  over  the  p 
future  inputs,  i .e. , 
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u(t)  =  u(t+l)  =  ...  -  u(t+p-l) 


Then  the  problem  is  reduced  to  the  determination  of  the  single  value 
u(t)  such  that  (4.21)  is  minimized.  It  is  clear  that: 


F  i  2 

min  l  w^y  (t+i+1)  -a  y(t)| 

(u(t),. . . »u( t+p-l) }  1-0 

p-1  .•  o 

<  min  l  w.|y  (t+i+1) -a  y(t)r 

{ u ( t )  =  ...  =  u( t+p-l) }  1-0 

Moreover,  in  the  minimization  shown  in  the  right-hand  side  of  the  above 
inequality,  the  optimum  value  of  u(t)  depends  on  the  length  p  of  the 
horizon  of  prediction,  while  in  the  left-hand  side  the  optimum  value  of 
the  input  sequence  fu(t),  u(t+l),,..,  u( t+p-l) }  is  independent  of  p. 
Schematically  the  differences  are  shown  in  Figure  19. 


Figure  19.  Effect  of  constant  inputs  on  optimization. 


l 
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Another  important  property  is  that  the  stability  of  the  system  is 
going  to  be  influenced  not  only  by  the  reference  model  (via  the  choice 
of  a)  but  also  by  the  p  parameters  w^.  It  will  be  seen  that  proper 
choice  of  w^.'s  (i  =  0,...,p-l)  can  stabilize  a  nonminimum  phase  system. 
Again,  it  is  worth  noting  that  the  introduction  of  the  horizon  p  and 
the  parameters  wi  (i  =0,...,p-l)  allowing  the  stabilization  is  done  at 
the  cost  of  degrading  the  distance  minimization  between  the  reference 
model  and  the  actual  output  during  a  transitional  period. 

In  the  following  two  subsections  we  treat  the  case  of  p  =  2.  The 
results  are  easily  extendable  to  p>2. 


4.6.1  Two  Step  Ahead  Prediction 

Consider  the  particular  case  where  p  =  2  and  W2=l.  That  is,  at 
period  t,  we  are  seeking  to  determine  the  inputs  u(t)  and  u(t+l)  such 
that  the  predicted  yp(t+2)  matches  perfectly  the  desired  reference  yr(t+2) 


yr(t+l)  =  ay(t)  (we  assume  that  c =  0) 
yr(t+2)  =  a2y(t) 


►  reference  model 


(4.22) 


y  (t+1)  =  y(t)  +  hT[u(t)  -  u(t-l)]l 
v  f  prediction  model 

yp(t+2)  =  y(t)  +  h  [u(t+l)  -  u(t)]J 


y(t+i)  =  hTu(t)  ) 

j  >  plant 
y(t+2)  =  h  u(t+l)J 


The  equality  condition  becomes: 


(4.23) 


(4.24) 
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(4.25) 


yr(t+2)  =  yp(t+2) 

That  is, 

h1u(t+l)  +  h2u(  t)  =  {-(l-a2)hT  +  hT(I-U2)}u(t-l) 

~  o 

where  hj,  h2  are  the  first  two  elements  of  h,  and  U  is  the  two  step 

ahead  shift  matrix.  The  additional  assumption  that  u(t)=u(t+l)  (constant 
input  over  the  prediction  horizon)  implies  that: 

u(t)  =  u(t+l)  =  ~7~  {-(l-a2)bT  +  hT(I-U2)}u(t-l)  (4.26) 

hl+h2 

N  can  be  assumed  even  (if  it  is  not,  we  can  always  consider  a  new 
impulse  response  of  length  N+l  with  hN+j  =  0) .  Therefore,  N  =  2q.  Then 
let  us  define  a  new  vector  v  of  dimension  N/2  such  that: 

v(t)  =  u(t)  =  u(t+l) 
v(t+l)  =  u(t+2)  =  u(t+3) 

v(t+q-l)  =  u(t+2q-2)  =  u(t+2q-l) 

Thus  equation  (4.25)  becomes: 

(hj  +  h2)v(t)  =  [hj  +  h2  -  (l-a2)(h1  +  h2)]v(t-l) 

+  Cho  +  R4-(l-a2)(h-  +  h4)]v(t-2) 

34  34  (4.27) 

+  -  .  . 

j  +  h^  -  ( 1  ~ci  )  (b^_j  +  h^)]v(t_q) 
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Define: 


r 

C\J 

-C 

+ 

rH 

-C 

\ _ 

-  +  fi2- 

h3  +  h4 

and  g  = 

h3th4 

z' 

XI 

+ 

rH 

1 

>-ltSN 

Then  the  above  equations  can  be  written  as: 


9xv(t)  =  [gT( I-U)  -  (l-a2)gT]v(t-l) 

=  (l-«2)9i3v(t-l)  +  ...  +  [gq-  (l-ot2)gq]v(t-q) 


(4.J 


where  ^  and  g^  are  components  of  g  and  g,  or: 

gTv(t)  =  ot2gTv(t-l)  +  gT[v(t)  -  y(t-l)]  (4.29) 

The  reader  familiar  with  MAC  formulation  will  easily  see  the  simi¬ 
larities  of  equations  (4.27)  and  (4.28)  with  equations  (4.8)  and  (4.10) 
describing  the  one-step  prediction  scheme.  Equations  (4.27)  and  (4.26) 

can  be  deduced  from  (4.8)  and  (4.10)  by  replacing  h,  h  and  a  with  g,  g 
2 

and  a  ,  respectively.  The  difference  is  in  the  dimension  of  the  system 
(N/2  instead  of  N).  All  the  considerations  of  Sections  4.4  and  4.5  can 
be  applied  to  the  present  system  of  dimension  N/2.  That  is,  if  gT  is 
minimum  phase,  then  one  would  choose  Then  equation  (4.26)  will 

yield  a  stable  input  sequence  v(t)  and  by  equations  (4.22)  and  (4.24),  we 
will  have 
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(4.30) 


y(t+2)  =  yr(t+2)  =  Ci2y(t) 


while 

y(t+i)  t  yr(t+l)  =  ay(t) 

Then  by  equation  (4.30), 

lim  y(t+2)  =  lim  yr(t+2)  =  0 

t"**5  t-X” 

and  by  equation  (4.23), 

1 im  y(t+l)  =  1 im  y(t)  +  hT  lim  [u(t)  -  u(t-l) J  =  1 im  y(t)  +  0 

t-*»  t-*®  t-*®  t-*°° 


Hence, 


lim  y(t+l)  =  lim  y(t)  =  lim  yr(t)  =  0 

t-KX>  t-*00  t-*00 

That  is,  although  during  a  transitional  phase  the  perfect  matching  of  yr 
and  y  occurs  each  two  steps  (t+2,  t+4,  ....  t+2j),  as  t  increases  y(t;  _»n- 
verges  to  a  desired  value  c  (here  c  =  0) .  Figure  20  shows  such  a  behavior 
for  a  6-order  system. 

Now  if  the  vector  corresponds  to  a  nonminimum  phase  system,  then 
one  has  to  determine  an  N/2-dimensional  vector  stabilizing  the  dynamics 
of  (4.28).  Both  methods  of  Sections  4.4  and  4.5  can  be  used.  The  clear 
advantage  is  that  the  dimension  N  has  been  reduced  by  half.  It  is  easy 
to  see  that  a  p-steps  ahead  prediction  scheme  would  reduce  the  dimension 
of  the  system  to  be  stabilized  to  N/p  instead  of  N. 
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1.,  -2.2,  2.18,  -.4420,  -1.4303,  .6594].  a=  .8,  c  =  0. 

Figure  20.  Two-step  ahead  prediction  for  the  nonminimum  phase. 


One  can  formulate  the  p-steps  ahead  prediction  problem  with  wi  =  0 
for  i  <  p  and  wp  =  1  as  a  one-step  ahead  prediction  problem  where  the  N- 
dimensional  previous  impulse  response  hT  is  represented  by  a  new  (N/p)- 
dimensional  impulse  response  g^. 

4.6.2  Weighted  Distance 

In  the  sequel  we  will  attempt  to  see  how  the  choice  of  the  weighting 
factor  w.  of  equation  (4.21)  can  be  used  to  stabilize  a  nonminimum  phase 
system.  Again,  for  the  sake  of  simplicity  we  consider  a  horizon  of  pre¬ 
diction  of  length  p  =  2.  We  seek  an  input  u(t),  constant  over  two  periods, 
which  minimizes  the  distance 

J2  =  wlCyp(t+1)  ‘  yr(t+1):]2  +  W t+2)  '  yr(t+2)]2  (4-31) 

where  y,  yr  and  yp  are  subjected  to  equations  (4.22)  through  (4.24).  We 
can  normalize,  without  loss  of  generality,  w^  and  by  requiring  that 

wl  +  w2  -  1 

The  input  u(t)  which  minimizes  Jg  is  the  solution  of  the  following  equa¬ 
tion,  where  u(t+l)  =  u(t): 

3 

3uft)  =  0 

which  leads  to: 
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0  =  u(t+l)[w1h12  +  w2(h1  +  h2)2]  +  fiT[w1h1(U-I)  +  (hj  +  h2)w2(U2-I)]u(t-l) 


+  y(t)[h1(l-a)w1 +  (hj  +  h2)(l-a2)w2] 


(4.32) 


or 

0  =  tKt+lKwjhj2  +  w2(hj  +  h2)2]  +  Wjhj[h^(U-I)  +  (l-a)h^]u(t-l) 

+  w2(hj +  h2)[hT(U2-I)  +  (l-a2)hT]u(t-l)  (4.33) 

The  above  equation  determines  the  input  u{t)  which  minimizes  J2,  as 
a  function  of  the  past  input  vectors  u(t+l).  It  is  easy  to  verify  that 
for  Wj  =  0,  equation  (4.33)  becomes  identical  to  equation  (4.27)  which 
corresponds  to  two-step  ahead  prediction.  Conversely,  for  w2 =  0  equa¬ 
tion  (4.33)  is  identical  with  the  equation  corresponding  to  the  one-step 
ahead  prediction  (see  Section  4.3). 

As  in  the  previous  sections,  the  controller  can  always  determine 
(N+l)  parameters  (wj,  hj,...,  h^)  such  that  the  characteristic  equation 
corresponding  to  equation  (4.33)  has  N-I  predetermined  poles.  However, 
this  solution  presents  no  advantage  over  the  ones  described  in  Section  4.4, 
and  moreover  it  does  not  specify  the  role  of  w^. 

Let  us  require  that  the  nathematical  model  represent  the  plant  exactly; 
that  is,  hT  =  hT.  Then  there  remains  only  one  parameter  Wj  to  adjust  'or 
the  stabilization  of  the  system.  By  equating  h^  to  hT  in  equation  (4.33), 
we  deduce: 
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it-- 


0  =  +#2(112  +  l^)^]  +  h^lWjhjCU-al]  +  #2  ( hj  +  ^KU^-a^IjJuft-l) 

(4.34) 

Clearly,  it  is  not  possible  to  bring  all  the  unstable  poles  of  the  above 
dynamic  to  predetermined  values  inside  the  unit  circle  by  adjustment  of 
only  one  parameter  v^.  However,  it  might  be  possible  to  stabilize  the 
system  with  a  proper  choice  of  W2;  that  is,  to  bring  the  unstable  poles 
inside  the  unit  circle.  It  is  clear  that  even  this  latter  attempt  is  not 
always  assured  of  success,  particularly  when  the  number  of  unstable  poles 
of  H(z)  increases.  However,  experience  has  shown  that  in  those  cases,  by 
increasing  sufficiently  the  length  p  of  the  horizon  of  prediction,  it  be¬ 
comes  possible  to  stabilize  the  system.  An  intuitive  explanation  of  the 
above  phenomenon  is  that  by  increasing  the  length  p,  one  introduces  p-1 
additional  parameters,  w^ ,  which  then  can  be  adjusted  to  stabilize  the 
system.  Comparing  the  above  observation  with  the  approach  of  Section  4.4, 
where  the  horizon  of  prediction  is  equal  to  1,  it  appears  that  the  choice 
of  hT  or  wi  are  complementary  in  the  sense  that  either  one  fixes  the  horizon 
of  prediction  at  p =  1  and  thus  has  to  determine  the  appropriate  vector 
h  ,  or  one  fixes  h  (h  = h  )  and  thus  has  to  determine  the  appropriate 
length  p  and  the  weights  w^  (i  =  l,...,p). 

The  practical  implementation  of  the  above  scheme  needs  minor  altera¬ 
tion  in  the  existing  MAC  software.  But  the  proper  choice  of  w.j  needs 
much  acquaintance  and  experience  with  the  particular  system  to  be  sta¬ 
bilized.  Figure  21  displays  an  example  of  such  stabilization  for  a  system 


of  order  5. 
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4.7  Summary 

The  problem  of  nonminimum  phase  systems  has  been  posed  within  the 
framework  of  MAC.  The  theoretical  analysis  has  shown  that  a  rather  wide 
range  of  solutions  for  the  control  of  such  systems  exists.  These  dif¬ 
ferent  types  of  solutions  are  closely  related  to  each  other.  Moreover, 
their  implementation  needs  minor  alterations  of  the  regular  MAC  software. 

In  practice,  the  specific  solution  to  choose  depends  on  the  particular 
system  and  also  on  the  desired  response  behavior,  and  generally  there  is 
no  one  "once-through"  design  method  yielding  an  implementable  solution. 

In  most  cases,  the  proper  solution  is  the  outcome  of  a  more  or  less 
lengthy  design  process  involving  both  theoretical  considerations  and  also 
simulation  results. 
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SECTION  V 


MATHEMATICAL  FORMULATION  AND  ANALYSIS  OF 
MAC  OF  CONTINUOUS  TIME  SYSTEMS 


Overview 

In  this  section  we  develop  a  mathematical  formulation  for  the  MAC 
scheme  applied  to  continuous  time  systems  rather  than  to  discrete  time 
systems  as  discussed  in  Section  III.  This  framework  provides  a  pre¬ 
cise  formulation  of  the  main  components  of  the  control  scheme  and  clearly 
indicates  the  many  various  possibilities  for  implementing  the  basic  con¬ 
trol  strategy.  In  addition,  this  framework  provides  a  means  of  analyzing 
the  behavior  of  the  MAC  scheme  in  different  situations,  and  permits  an 
assessment  of  its  computational  requirements.  The  section  is  divided 
into  three  sections:  the  first  section  presents  the  general  formulation; 
the  second  discusses  the  optimal  solution  and  the  form  of  the  MAC  control 
law;  and  the  third  section  briefly  discusses  the  computational  aspects 
of  implementing  the  MAC  control  strategy. 

5.1  General  Formulation  of  MAC  Strategy 

Simply  stated,  the  MAC  philosophy  of  control  is  to  choose  fictitious 
future  controls  in  order  to  minimize  a  weighted  quadratic  distance  be¬ 
tween  a  predicted  output,  which  is  a  function  of  the  fictitious  future 
controls,  and  a  desired  future  output  (reference  trajectory).  The  first 
step  or  first  few  steps  of  the  fictitious  future  control  are  actually 


applied,  and  after  a  short  interval  of  time  a  new  fictitious  future  con¬ 
trol  is  computed  and  the  process  repeats  itself.  To  formulate  this  pre¬ 
cisely,  let  us  make  the  following  definitions. 

Predicted  future  output:  yp(s,t)  is  the  output  predicted  at  future 
time  t+s  when  t  is  the  present  time.  Note  that  s>0.  The  predicted  future 
output  yp(s,t)  depends  on  actual  outputs  and  inputs  up  to  the  present 
time  t  and  on  the  fictitious  future  inputs  up  to  the  future  time  s+t. 

Fictitious  future  inputs:  u^(s,t)  denotes  the  fictitious  future 
input  supposed  to  be  applied  at  the  future  time  t+s  when  t  is  the  present 
time. 

Desired  future  output  ( reference  trajectory):  yr(s,t)  denotes  the 
output  desired  at  the  future  time  t+s  when  t  is  the  present  time.  This 
reference  trajectory  depends  on  actual  outputs  and  inputs  up  to  the  pre¬ 
sent  time  t.  It  may  also  depend  on  a  prespecified  objective  function  or 
desired  set-point. 

The  MAC  strategy  is  to  minimize  a  weighted  quadratic  distance  between 
yp(*,t)  and  yr(*,t)  of  the  form 

J(uf(-,t))  =  |yp(s,t)  -  yr(s,t)|2w(s)  ds  (5.1)6 

where  J  depends  on  u^(*,t)  through  yp(*,t).  In  minimizing  J  we  restrict 
uf(*»t)  to  a  fixed  set  of  admissible  future  inputs.  The  actual  input 
u(t)  applied  at  time  t  is  simply 

^In  this  discussion  the  input  and  output  values  are  finite  dimensional 
vectors.  The  norm  |v|  is  the  usual  Euclidean  norm  for  the  vector,  i.e., 

2  2  2 

j v |  =v  j  +  .  ..  +  v^j  where  v  is  d-dimensional  and  are  the  coordinates 
of  v  with  respect  to  fixed  coordinate  basis. 
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u(t)  =  uf*(0,t) 


(5.2) 


where  uf*{',t)eUf  minimizes  J. 

To  illustrate  the  control  strategy  more  clearly  we  consider  typical 
implementations.  First,  we  need  to  assume  a  particular  functional  re¬ 
lationship  between  the  fictitious  future  control  and  the  predicted  future 
output.  In  other  words,  we  need  a  model  for  the  future  response  of  the 
plant.  If  we  choose  a  linear  model  (as  we  usually  do),  then  yp  can  be 
expressed  conveniently  as  the  sum  of  the  zero  input  response  ypzi  and 
the  zero  state  response  ypzs  of  the  system 

yp(s,t)  =  yp2i(s,t)  +  ypzs(s,t)  (5.3) 

That  is,  y  .(s,t)  =  yp(s,t)  if  the  future  input  is  zero  (uf(s,t)=0  for 
all  s>0).  The  zero  state  response  ypzs  then  gives  the  future  output  due 
to  future  input  alone. 

We  assume  that  ypz$  is  given  by  an  impulse  response  model  in  terms 
of  u^  as  fol lows: 

y_7<.(s,t)=f  H(s-o,t)  uf(o,t)  do  (5.4) 

pzs  Jq 

In  the  following  discussion  we  will  also  assume  that  H(s-o,t)  = H(s-o)  is 
independent  of  the  current  time  t. 

There  are  several  possibilities  for  representing  the  zero  input 
response  y  . .  For  example,  one  can  use  the  impulse  response  of  the 
actual  past  inputs: 
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(5.5) 


ypzi ( 5 *t)  =  {  H(s-o)u(ct)  da 

Alternatively,  one  might  use  a  state  model  of  the  form: 

ypzi(s,t)  =  F(s)  *(L)  (5-6) 

where  x(t)  is  a  function  of  actual  past  inputs  and  outputs  computed  by 
some  finite  dimensional  filter.  Other  possibilities  (e.g.,  ARMA ‘type 
models)  exist,  and  one  can  also  consider  combinations  of  these.  The 
important  thing  to  note  is  that  the  computations  of  ypzi  and  ypz$  can 

be  carried  out  independently.  On  the  one  hand,  it  seems  most  computa¬ 
tionally  efficient  to  use  the  impulse  response  model  (5.4)  to  compute  ypzs- 
On  the  other  hand,  the  impulse  response  model  may  not  be  the  most  ef¬ 
ficient  procedure  for  computing  y  .,  and  it  is  a  good  idea  and  easy  to 
tailor  the  zero  input  response  model  to  the  particular  problem. 

There  are  also  several  choices  for  the  reference  trajectory  y  .  Un¬ 
fortunately  (for  purposes  of  mathematical  analysis),  these  choices  are 
made  more  on  the  basis  of  intuition  than  systematic  theory.  However,  it 
is  easy  to  express  the  most  conron  choices  in  mathematical  terms,  and  we 
can  subsequently  analyze  some  of  the  implications  of  these  choices.  The 
most  common  choice  is  the  first-order  reference  trajectory  of  the  form 

yr(s,t)  =  e"s/Ay(t)  +  (1  -  e"s/A)c(t)  (5.7) 

where  A  is  a  specified  time  of  response,  y(t)  is  the  actual  input  at 
time  t  and  c(t)  is  a  set-point  or  objective  at  time  t  (and  c(t)  is  often 
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constant  in  this  case).  Another  choice  of  the  reference  trajectory  is 
given  by 


yr(s,t)  =  c(t+s)  (5.8) 

where  c(*)  is  a  pre-specified  objective  function.  (Here  c(t)  usually 
tends  to  a  constant  value  as  t-*-°°.)  This  choice  corresponds  to  the  out¬ 
put  predictive  schemes  of  Reid  et  al_.  (1979,  1980). 

Let  yQ(s,t)  denote  the  difference 

yQ(s,t)  =  yr(s,t) -ypzi(s,t)  (5.9) 

At  time  t,  we  must  minimize 


J(uf(*,t)) 


J0  lyo(s’t}  - 


H(s-a)uf(o,t)da|  w(s)  ds 
0  T 


(5.10) 


with  respect  to  u^(*,t)eUf.  The  set  of  future  admissible  controls 
must  be  a  subset  of  a  finite  dimensional  vector  space  of  inputs  (since 
we  have  to  do  the  optimization  numerically  on-line)  and  it  must  be  such 
that  J  has  a  unique  minimum  over  U^.  Here  we  assume  that  U ^  can  be  ex¬ 
pressed  as  a  set  of  unique  linear  combinations  of  input  functions  un(s) 
defined  for  s>0. 


uf(*,t)  =  l  v.  uV)  (5.11) 

T  i  =  l  1 

where  v =  (v^.Vg, . . . ,vn)^ e  V  and  V  is  a  closed,  bounded  polyhedron  in 
Fn.^  The  optimization  of  J  is  performed  inFn  on  the  finite  dimensior'l 
^Here  and  elsewhere  (  )  denotes  the  transpose  of  both  vectors  and  matrices. 


vectors  V.  In  terms  of  veV,  J  is  a  quadratic  form  in  v. 


J(v)  =11  P1kVk  -  2  l  q  (t)v. 
j=l  k=l  JK  J  K  j=l  3  J 


(5.12) 


where 


qj(t)  =  |°°  w(s)yQ(s,t)T|  HCs-aJu^CoJdaJds 


(5.13) 


and 


u^(a)^H(s-o)^H(s-T)u*c(x)dodrlds 


(5.14) 


where  j,k  =  1,2,. ..,n.  Note  that  Pjk  is  independent  of  current  time  t 
and  doesn't  need  to  be  computed  on-line.  In  order  to  ensure  a  unique 
minimum,  we  assume  u’(*)  are  chosen  so  that  [pjk]  is  a  nonsingular  matrix 


(hence  positive  definite). 

The  linear  terms  qj(t)  do  depend  on  t  through  yg(s,t)~that  is,  due 
to  the  predicted  zero  input  response  y  .(»,t)  and  the  reference  trajec¬ 
tory  y  (*,t).  However,  we  can  often  arrange  the  computation  of  q(t)  so 

■  J 

that  most  of  the  calculations  are  done  off-line.  For  example,  if  ypZ_.  is 
given  by  (5.6)  and  y  is  given  by  (5.7),  then  q.(t)  can  be  written 

'  J 

qj(t)  .  *  y(t)T„/  +  c(t)Tqj3  (5.15) 


where  q^ 

J 


1 


are  vectors  which  can  be  computed  off-line  as 
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(5.16) 


Qj1  =  -j^  w(s)F(s)Tjj  H(s-o)u^(a)do|ds 
=  w(s)  e  S/^jj  H(s-o)uJ (o)dojds 

=  w(s)  (1-e  S/,A)||  H(s-a)uJ (a)dajds 


(5.17) 


(5.18) 


At  this  point  let  us  review  the  main  steps  of  the  control  algorithm: 


MAC  Strategy 

(i)  Compute  certain  terms,  such  as  P..  in  (5.14)  and  q1?'  in 

J  K  J 

(5.16)  -  (5.18)  off-line. 

(ii)  At  each  control  update  time  t,  compute  q.(t)  from  (5.13) 
or  (5.15). 

(iii)  Minimize  the  quadratic  form  J(v)  in  (5.12)  with  respect  to  v 
and  its  linear  inequality  constraints.  Let  v*(t)  denote 
the  minimum  solution. 

(iv)  Apply  the  actual  control  input 


(5.19] 


u(t)  =  l  v.*(t)uK(0) 
k=l  K 

up  until  the  next  control  update  time . ^ 
(v)  Return  to  (ii)  and  repeat. 


5.2  Solution  of  Optimization  and  Form  of  MAC  Control  Law 


The  type  of  control  law  generated  by  the  MAC  strategy  depends  on  the 
way  in  which  the  terms  vk*(t)  in  (5.19)  depend  on  the  past  data.  In  turn, 

Q  ^  L 

More  generally,  one  would  have  u(t+s)=  l  v£(t)u  (s)  for  times  t+s  up 
until  the  next  update.  k=l  K 


the  dependence  of  Vj,*(t)  on  the  past  data  is  determined  by  the  dependence 
of  q . ( t )  in  (5.12)  on  the  past  data.  Note  that  q.(t)  depends  linearly  on 

J  J 

the  past  in  the  cases  we  considered  in  the  previous  section.  If  the 
minimum  {v^*(t)}  is  unconstrained,  then  it  too  depends  linearly  on  the 
past.  Let  v* ( t )  denote  the  n-vector  with  coordinates  qk( t.)  and  let  P 
denote  the  n*n  matrix  [Pkj].  Then,  the  unconstrained  minimum  of  J(v) 
in  (5.12)  is  given  by 

v*(t)  =  P~*q(t)  (5.20) 

The  control  law  (5.19)  can  be  written 

u(t)  =  Dv*(t)  (5.21) 

u  th 

where  u  (0)  is  the  k~  column  of  the  matrix  U.  In  case  of  unconstrained 
v*(t)  we  get  the  linear  control  law 

u ( t )  =  UP_1q(t)  (5.22) 

If  the  reference  trajectory  is  given  by  (5.7)  and  the  zero  input  response 
is  given  by  (5.6),  then  (5.15)  is  true  and  (5.22)  can  be  rewritten  expli¬ 
citly  in  terms  of  the  past  data  as 

u(t)  =  UP'1(Q1x(t)  +  Q2y(t)  +  Q3c(t))  (5.23) 

o  p  T  th 

In  this  expression  Q  is  the  matrix  with  (q.  )  as  its  j-  row. 

In  case  v*(t)  is  constrained,  the  control  law  is  nonlinear  and  much 
more  complicated.  However,  something  more  can  be  said  about  the 
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functional  form  of  the  control  law.  Consider  the  figure  below  (Figure  22), 


In  this  case  n  =  2  and  the  constraints  on  v1  and  have  the  form 


vk  *  vk  *  vk  » 


k  =  1,2 


(5.24) 


The  constraint  set  V  is  the  square  region  indicated  in  the  figure.  The 
dashed  lines  in  the  figure  represent  directions  orthogonal  to  the  con¬ 
straints  v(<=vk±.  The  orthogonality  of  the  directions  is  with  respect 
to  P  and  in  general  these  directions  would  be  oblique  to  the  constraints. 
For  simplicity  we  have  drawn  them  perpendicular  (this  corresponds  to  a 
diagonal  P). 


Let  v**(t)  denote  the  unconstrained  minimum 


v**(t)  =  P_1q(t) 


(5.25) 


which  may  give  v**(t)  ^V.  It  is  not  hard  to  see  that  in  each  of  the 
regions  0,  1,  2,  3,  4,  5,  6,  7,  8  shown,  the  constrained  minimum  v*(t) 
is  an  affine  function  of  the  unconstrained  minimum  v**(t).  That  is  to 


say,  v*(t)  can  be  written  as 


v*(t)  =  $'  (v**(t)  )v**(t)  +  4>'  (v**(t) ) 


(5.26) 


where  <J>'(v**)  is  an  n  *  n  matrix  which  is  piecewise  constant  on  the  regions 
0-8,  and  <j»'(v**)  is  a  piecewise  constant  n-vector  on  the  same  regions. 
Thus,  we  say  that  the  general  MAC  control  law  given  by  (5.21),  (5.25)  and 
(5.26)  is  piecewise  linear-affine.  The  pieces  are  especially  nice,  being 
defined  by  simple  linear  inequalities.  Unfortunately,  as  n  grows,  the 
,number  of  pieces  grows  exponentially.  Nevertheless,  there  may  be  effective 
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methods  of  anaTyzing  this  type  of  control  law.  For  instance,  it  is  an 
example  of  the  variable  structure  systems  studied  by  Utkin  (1977).  Also 
in  certain  cases  it  might  be  possible  to  analyze  the  behavior  of  v*(t) 
effectively  in  terms  of  the  unconstrained  minimum  v**(t).  However,  it 
must  be  remembered  that  in  general  v**(t)  will  evolve  according  to  a 
nonlinear  law  also.  When  the  inequality  constraints  are  not  as  simple  as 
in  (5.24)  (for  example,  if  there  are  state  constraints),  the  situation 
becomes  much  more  complex  than  pictured  in  Figure  22.  It  is  still  possible 
to  express  the  control  law  as  piecewise  linear-affine,  but  the  pieces  be¬ 
come  much  more  complicated.  In  any  case,  we  will  not  pursue  the  inves¬ 
tigation  of  the  nonlinear  control  law  any  further  here.  Appendix  A  inves¬ 
tigates  the  nonlinear  type  control  law  further  for  single  input-single 
output  systems. 

So  far  we  have  seen  that  the  general  MAC  control  law  is  nonlinear  and 

q 

quite  difficult  to  analyze.  The  MAC  control  law  is  linear  in  the  case 
that  there  are  no  inequality  constraints  on  the  set  of  admissible  future 
controls  so  that  V=Fn.  However,  the  linear  case  offers  its  own  dif¬ 
ficulties  due  to  the  possibility  of  unbounded  control  inputs  (since  J 
does  not  weight  controls  at  all).*®  The  choice  of  and  the  reference 
trajectory  in  this  case  can  make  the  difference  between  stable  and  un¬ 
stable,  robust  and  non-robust  control. 

Before  finishing  this  section,  we  note  that  the  most  common  choice 

of  Uf  is  the  set  of  inputs  piecewise  constant  over  disjoint  intervals 

of  time.  Thus,  the  basis  inputs  un(*)  have  the  form 

g 

although  the  particular  piecewise  linear-affine  structure  offers  some 
promise  of  analysis. 

10Control  weights  are  easily  added  and  easily  analyzed,  of  course,  but 
their  use  has  not  been  required  in  the  examples  we  have  treated. 
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U1  (s)  =  U1,  S6l. 

=  0,  s*I. 


(5.27) 


where  u1  is  a  constant  vector  and  L  is  a  bounded  interval  in  [0,®).  The 
choice  of  1^  (the  so-called  "time  blocks"  or  "blocking")  greatly  affects 
the  type  of  unconstrained  control  computed  by  the  MAC  control  law.  Such 
a  choice  of  admissible  future  controls  leads  one  naturally  to  analyze 
the  problem  in  terms  of  discrete  time  models  as  presented  in  Section  III. 


5.3  Computational  Considerations 

In  case  of  input  or  output  constraints,  one  must  use  a  numerical 
optimization  algorithm  to  minimize  the  quadratic  form  (5.12)  with  respect 
to  the  constraints.  Currently,  a  simple  gradient  projection  algorithm 
(Canon,  Callum  and  Polak,  1970)  is  used  to  solve  this  optimization  problem. 
The  problem  is  set  up  as  a  quadratic  optimization  with  linear  inequality 
constraints  on  the  vector  v  appearing  in  (5.12).  Note  that  if  there  are 
constraints  specified  on  the  input  u(*)*  then  these  usually  show  up  as 
fairly  simple  constraints  on  v.  Output**  constraints,  on  the  other  hand, 
show  up  as  more  complicated  linear  inequality  constraints  on  v,  and  these 
tend  to  slow  down  the  computation  of  the  optimal  v*.  If  only  input  con¬ 
straints  are  present,  the  complexity  of  the  optimization  problem  depends 
only  on  the  dimension  n  of  the  matrix  P  in  (5.12).  In  most  cases  this 
dimension  is  small  (n<20).  However,  output  constraints  can  seriously 
increase  the  complexity  of  the  problem  even  if  P  has  small  dimension.  In 
this  case  it  would  be  advisable  to  use  a  small  number  of  constraints  on  v 

**State  constraints  can  be  formulated  as  (unobserved)  output  constraints 
in  this  framework. 
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to  approximate  the  desired  output  constraints. 

Another  important  aspect  of  computation  is  how  much  must  be  computed 
on-line  and  updated  each  time  one  wants  a  new  control.  As  noted  above, 
certain  quantities  such  as  and  q*j  in  (5.14),  (5.16)  -  (5.18)  can  be 
computed  off-line  and  stored  for  on-line  use.  Input  constraints  appear 
as  simple  linear  inequalities  on  the  v.  which  are  usually  independent  of 

J 

time.  Output  constraints,  on  the  other  hand,  may  be  time-dependent  when 
expressed  in  terms  of  linear  inequalities  in  the  variables  v.. 

J 

Although  there  are  a  number  of  sophisticated  algorithms  for  solving 
quadratic  optimization  problems  with  linear  inequality  constraints,  we 
are  more  concerned  at  present  with  obtaining  a  clear  understanding  of 
the  behavior  of  MAC  controlled  systems.  For  this  reason,  we  have  re¬ 
stricted  the  current  research  to  a  simple  gradient  projection  algorithm. 
Nevertheless,  it  is  good  to  note  that  this  simple  algorithm  is  sufficiently 
fast  to  update  the  control  inputs  for  the  problems  considered  so  far. 
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SECTION  VI 


MAC  FOR  CONTINUOUS-TIME  SYSTEMS  WITH 
DISCRETE  OBSERVATIONS  AND  UNOBSERVED  OUTPUTS 


Overview 

In  this  section,  we  investigate  MAC  for  an  incompletely  observed 
continuous  linear  time  invariant  system  containing  both  process  and  measure¬ 
ment  noise.  The  performance  index  is  taken  to  be  the  expected  value  of 
a  quadratic  function  of  the  continuous  tracking  error  and  control  effort. 

In  accordance  with  LQG  theory,  the  subject  is  analyzed  as  two  essentially 
separable  problems:  the  estimation  and  the  control.  The  estimation, 
which  is  of  minimum  variance  type,  is  done  by  Kalman  filtering.  For  the 
control  part,  two  types  of  solutions  are  considered.  One  is  based  on  the 
usual  MAC  approach  utilizing  impulse  response  functions  and  the  other  is 
based  on  state  feedback  control.  The  latter  approach  is  taken  here  to 
study  the  sample  rate  problem,  i.e.,  determination  of  the  effect  of  sample 
rate  on  controller  performance.  The  organization  of  this  section  is  as 
follows.  The  problem  is  formally  stated  and  discussed  in  Section  6.1. 
Section  6.2  analyzes  the  optimal  estimation  problem.  Given  the  output  of 
the  estimator,  the  optimal  control  is  discussed  in  Section  6.3.  Section  6.4 
is  devoted  to  analyzing  the  effect  of  sample  rate.  The  plant  model  in 
this  section  Is  assumed  to  be  in  a  state  vector  form.  This  is  convenient 
for  studying  the  problem  of  unobserved  outputs  and  sample  rate  selection. 

In  addition,  for  aerospace  vehicles,  state  vector  models  are  generally 


available.  MAC  implementation,  however,  is  not  restricted  to  purely 
impulse  response  models  or  state  vector  models.  As  shown  in  Section  V, 
it  is  possible  to  use  the  state  model  for  estimation-prediction  purposes 
and  the  impulse  response  model  for  control  computations.  Since  these 
approaches  are  mathematically  equivalent  for  the  idealized  case  of  un¬ 
constrained  inputs,  we  will  formulate  the  sample  rate  selection  problem 
in  terms  of  the  state  vector  approach,  thereby  exploiting  the  results  of 
LQG  theory. 


6.1  Statement  of  the  Problem 

Consider  a  discretely  observed  continuous-time  linear  time-invariant 
system,  given  by 


x(t)  =  Fx(t) +  Gu(t)  + w(t) 
y(t)  =  Hx(t)  +  v(t) 


(6.1) 


where 


x(*)eRn,  yHeRm,  u(«)eRr; 

E[x(t0)]  £  x(tQ); 

E[(x(t0)-x(t0))(x(t0)-x(t0))T]  6  Rx(t0); 

F  e  Rn*n,  GeRnxr,  HeRmxn; 
w(*)eRn.  v ( - )  e  Rm ; 
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w(*)  and  v{*)  are  zero-mean,  uncorrelated  white  noise  vectors  with: 


E[w(t)wT(t+T)]  ^  Q( t) 6(t ) ; 

E [ v ( t ) vT ( t+T ) ]  £  R( t) 6( x ) ; 

Q(*)eRnxn,  Q(*)>0;  R(-)eRw,  R(-)>0. 

The  system  is  observed  at  discrete  times  t..,  i  = 

y(t*)  =  Hx(tj)  +  v(tt)  (6.2) 

It  is  desired  to  determine  the  optimum  control  vector  u(-)  which 
minimizes  the  cost 


J  =  Eth 


t()+T(  ll^(t)  -  *  (t)  ||2  +  p||u(t) |2)dt] 
t0 


(6.3) 


where  yr(*)  is  the  reference  trajectory.  It  is  assumed  that  yr(‘)  has 
first  order  dynamics,  and  that  the  setpoint  is  at  the  origin.  Then 


ir(t)  * 

and  Yr(t0)  «  y(tQ) 


(6.4] 


and  the  reference  trajectory  is  given  as  a  function  of  the  initial  obser¬ 
vation  and  time  as 


¥r(t)  =  *(t0)e 


-( ( 1/t) (t-tQ) ) 


(6.5] 


no 


or 


Yr(t)  = 

where 

-(t-tJ/T 

a  =  e 

The  positive  constant  p  in  equation  (6.3)  is  the  cost  associated  with  the 
input  energy. 

It  is  assumed  that  at  time  te[tg,tg+T],  the  controller  has  the  in¬ 
formation  Y(t)  =  {y(t|);  tg<t^t},  i.e.,  all  samples  prior  to  t.  The 
determination  of  u(t)  is  based  on  the  information  set  Y(t). 

It  is  well  known  that  the  determination  of  the  optimal  input  u ( • ) 
involves  an  estimation  problem  and  a  control  problem.  In  the  remaining 
part  of  the  present  section,  the  separation  property  (separation  theorem) 
is  derived  for  the  system  given  by  equations  (6.1)  and  (6.2). 

The  cost  J  can  be  rewritten  as: 

J  =  E{E[>2  Jto+T(  l!y(t)  -  yr(t)  | ^  +  p II u ( t )  ||2)dtj  Y(t0+T)]}  (6.6) 

0 

where  E[...  Y(t)]  denotes  the  conditional  expectation  given  an  information 
set  Y(  •) »  while  E{...}  denotes  the  expectation  over  all  possible  informa¬ 
tion  sets.  Then 

tn+T  . 

J  =  E{4  (E[|*(t)-*r(t)|f  |Y(t0+T)]+pE[||u(t)||Z(Y(t0+T)])dt}  (6.7) 

Jt0 
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The  conditional  mean  and  covariance  of  the  random  vector  y(t)  are 


given  by: 


E[y(t)|Y(t)]  ^y(t) 

E[(y(t)-y(t))(y(t)-y(t))T|Y(t)]  $  Py(t) 


E[|y(t)-yr  (t)||2|Y(t)]  =  ll£(t)-yr(t)||2  +  tr  cov[y(t)-v  (t)] 


(6.8) 


where  yr(*)  is  defined  similarly. 

The  control  being  independent  of  the  observation  noise: 


E[|u(t)l2|Y(t)]  =  ||u(t)||2 


(6.9) 


it  remains  to  evaluate  cov[y(t)  -  y  (t)]  using 


cov[y(t)  -yr(t)]  =  cov[y(t)]  +  cov[yr(t)]  -  2E[(y(t)  - y(t) )T(yr(t)  -yr(t))] 

(6.10) 


It  follows  that: 

J  =  E[J  °  (I2(t) -yjt)||2  +  p||u(t)!2)dt  +  J5trf  cov[y(t)  - y  (t)]dt] 
t0  0 

(6.11) 


where 
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T  T  ~2(  1-Ot)  (  t-tn) 

cov[y  (t)]  =  P  (t)  =  E[HTxT(0)x(0)H]e  0 

yr 


Py  (t)  =  HTPx(t0)He 


2(l-a)(t-t0) 


P  (t)  is  independent  of  the  estimation  and  control  strategy.  The  task 
yr 

hence  reduces  to  minimizing 


J  =  E[% 


( Il2(t)  -  yr(t)  l2  +  pflu(t)  I2  +  tr  Py(t)  -  2E[(y(t)  -  y(t)  )T 


'(^r(t)-ir(t))|Y(t)])dt] 


(6.12) 


where  P  (t)  denotes  the  covariance  of  y(t)  conditional  on  the  information 
Y(t). 

The  above  development  shows  that  two  problems  must  be  solved;  first, 
it  is  necessary  to  design  an  estimator  to  obtain  y(t): 

2(t)  =  E[*(t)|Y(t)] 

This  estimator  is  of  the  minimum  variance  type.  The  second,  in  a  sense 
independent,  problem  is  to  derive  the  control  u  that  minimizes  equation 
(6.12).  The  following  section  deals  with  the  design  of  the  estimator. 


6.2  Optimal  Estimation 

The  optimal  (in  the  minimum- variance  sense)  estimate  y(t)  of  the 
system  (6.1)  given  the  discrete  observations  y(tt),  i  =  l,...,p,  is  generated 


(see  Jazwinski,  1970  and  Bryson  and  Ho,  1975)  by  the  linear  filter: 
x(t)  =  Fx(t)  +  Gu(t)  +  K(ti)(y(tt)  -  Hx(tT))5(t-ti) 

(6.13) 

2(t)  =  Hx(t) 

for  te[t^,t^+^),  i  =  l,...,p.  Here  K( t^ )  is  the  Kalman  gain  and  the  term 
(y(tj)-y(tT))  =  (y(tt)-Hx(tT))  (6.14) 

is  the  innovation,  which  contains  the  new  information  due  to  the  observa¬ 
tions  at  time  t^ ,  ^(t*).  It  may  be  noted  that  in  general 

y(tT)  f  £(tT)  =  lim  Hx( t)  (6.15) 

t-*tT 

The  system  (6.13)  describes  the  evolution  of  the  estimate  y(t)  be¬ 
tween  instantaneous  observations  at  times  t.. .  The  Kalman  gain  K( t^ )  is 
obtained  in  terms  of  the  solution  of  a  Riccati  equation.  Since  these 
results  are  well  documented  elsewhere  (see,  e.g.,  Astrom,  1970,  Kwakernaak 
and  Sivan,  1972),  we  omit  the  details  here. 

6.3  Optimal  Control 

Once  the  complete  state  of  the  system  has  been  estimated,  all  the 
controlled  outputs  (observed  and  unobserved)  can  be  predicted  for  dif¬ 
ferent  fictitious  inputs  in  the  future  (see  Section  V).  One  may  use 
either  the  state  vector  model  or  the  impulse  response  model  for  predic¬ 
ting  the  zero  state  response,  depending  on  the  dimensionality  of  each 
representation  and  the  associated  computational  effort.  The  optimization 
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problem,  however,  is  more  efficiently  solved  using  a  discretized  impulse 
response  representation,  particularly  when  hard  constraints  have  to  be 
satisfied. 

The  computation  of  the  optimal  controls  using  the  impulse  response 
model  is  similar  to  that  carried  out  in  Section  5.1.  The  certainty 
equivalence  principle  (Astrom,  1970),  however,  does  not  apply  strictly 
in  the  constrained  control  (or  state)  case.  The  analysis  of  this  dif¬ 
ficult  case  requires  further  investigation.  We  assume  here  that  the 
effect  of  constraints  is  taken  into  account  indirectly  by  the  control 
effort  weighting  term  in  the  performance  index.  The  optimal  control 
problem  is  thus  equivalent  to  an  L.QG  problem,  which  may  be  stated  in 
conti nuous-time  or  discrete-time  depending  on  the  nature  of  the  control 
input.  Most  commonly,  the  control  input  is  kept  constant  (or  linearly 
interpolated)  between  sample  points  so  that  a  discrete-time  formulation 
is  appropriate.  Notice  that  the  sample  rate  for  the  control  input 
application  need  not  be  the  same  as  the  sample  rate  for  output  observa¬ 
tion.  The  former  restricts  the  class  of  inputs  over  which  optimization 
is  performed,  whereas  the  latter  influences  the  prediction  accuracy. 

The  effect  of  these  two  sample  rates  on  controller  performance  can  be 
studied  in  terms  of  the  solution  to  two  discrete-time  Riccati  equations, 
the  control  Riccati  equation  obtained  by  discretizing  the  system  equations 
and  the  performance  index  at  the  controller  sample  rate,  and  the  estima¬ 
tion  Riccati  equation  obtained  by  discretizing  the  system  equations  at 
the  observation  sample  rate.  The  relevant  equations  are  given  in 
Kwakernaak  and  Sivan  (1972),  pp.  542-550.  (See,  in  particular,  equations 
6-486,  6-487,  6-488,  6-523,  6-524  and  6-525.) 
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6.4  Summary 


It  was  shown  that  the  problem  of  MAC  design  for  stochastic  continuous¬ 
time  systems  with  noisy  incomplete  discrete-time  observations  can  be 
reduced  to  the  solution  of  an  estimation  problem  and  a  deterministic 
control  problem.  The  former  is  solved  in  terms  of  a  Kalman  filter  (see 
Section  VII  for  an  alternative  approach)  and  the  latter  is  shown  to  be 
solvable  either  using  LQG  theory  or  the  procedure  outlined  in  Section  V. 

A  study  of  the  sample  rate  selection  problem  is  indicated  in  terms  of  the 
solution  of  two  Riccati  equations. 
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SECTION  VII 


RECURSIVE  PREDICTORS 


Overview 

This  section  describes  the  preliminary  development  of  an  estimator 
to  work  with  IDCOM  in  the  control  of  systems  in  noisy  environments.  The 
results  of  Richalet  et  a1_.  (1978),  Mehra  et  al_.  (1978)  and  our  prelimi¬ 
nary  simulation  results  indicate  that  IDCOM  works  well  in  low-noise 
environments.  In  several  interesting  processes,  however,  this  ideal 
environment  cannot  be  guaranteed.  A  state  estimator  (Kalman  filter)  can 
of  course  be  included  in  IDCOM,  for  output  prediction,  when  a  good  state 
model  is  available.  But  the  requirement  of  a  good  state  model  is  a  severe 
one,  and  is  not  shared  by  the  control  optimization  part  of  IDCOM.  We 
were  thus  led  to  consider  alternate  predictors  which  are  compatible  with 
IDCOM  yet  offer  good  performance  when  noise  and  model  uncertainty  are 
present. 

This  section  begins  with  a  description  of  the  current  predictor 
in  IDCOM  (Section  7.1).  Section  7.2  then  discusses  a  simple,  useful 
modification  which  improves  performance,  and  Section  7.3  compares  this 
predictor  with  the  usual  one.  Autoregressive  Moving  Average  (ARMA)  models 
are  briefly  discussed  in  Section  7.4,  and  optimal  estimation  is  related 
to  this  formulation  in  Section  7.5. 


7.1  I DCOM  Prediction 

It  is  useful  to  note  how  IDCOM  currently  computes  the  state  (output) 
predictions  that  it  requires  for  control  computation.  The  IDCOM  fun- 
tions  may  be  outlined  as  in  Figure  23.  The  current  output,  y  ft),  i>  fed 

a 

to  two  IDCOM  blocks,  one  of  which  calculates  the  reference  trajectory 
for  several  future  time  steps  y^s.t)12  and  one  which  calculates  the  pre¬ 
dicted  output  for  zero-input  in  the  future  (i.e.,  no  future  controls 
applied)  ypzi(s,t). 


I  Desired 


Impulse 


Figure  23.  IDCOM  functional  diagram. 

Tp  ' 

This  notation  is  the  same  as  that  of  Section  V. 
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The  difference  between  these  terms  is  the  sequence  e(s,t)  of  future 
errors  to  be  removed  by  the  future  controls,  u(r,t).  The  only  control 
actually  applied  is  u(0,t),  the  current  choice,  since  the  calculation 
is  repeated  after  the  next  step. 

The  internal  (to  I DCOM )  impulse  response  model  is  usually  used  in 

two  places--the  control  calculation  and  the  zero-input  prediction.  These 

two  functions  do  not  require  the  same  model,  and  in  many  cases  an  alternate 

prediction  model  (e.g.,  for  noisy  environments)  or  control  model  (e.g., 

for  nonminimum  phase  systems)  may  be  desirable. 

In  general,  however,  IDCOM  uses  the  same  impulse  response  model, 

identified  (off-  or  on-line)  for  both  functions.  In  order  to  avoid  the 

dangers  of  "open  loop"  prediction,  however,  where  only  the  past  inputs 

and  internal  model  are  used,  IDCOM  modifies  the  standard  prediction  by 

adding  a  bias-compensation  term  to  improve  performance.  Thus,  where  the 

13 

"open  loop"  prediction  for  n-steps  ahead  would  be 


yozi  (n+1’t) 
pz10L 


T-n-1 

J0 


the  "closed  loop"  (bias  compensated)  version  is 


(n+l,t) 


(7.1) 


where 


y 


T 

j  HjU(t-i ) 


(7.2) 


13 


This  assumes  that  only  T  impulse  response  terms  are  stored. 
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1 


is  the  estimated  current  output  based  on  the  past  inputs  and  the  internal 
model.  Intuitively,  the  "closed  loop"  predictor  applies  the  bias  between 
the  actual  and  currently  predicted  outputs  to  all  future  outputs.  This 
technique  has  worked  quite  well  for  standard,  stable  plants  (with  finite 
impulse  response). 

7.2  Modification  for  Unstable  (Integrator)  Plants 

While  investigating  the  application  of  IDCOM  to  simple  systems  (e.g., 

integrators)  with  infinite  impulse  responses,  it  was  discovered  that  a 

slight  modification  to  the  above  scheme  resulted  in  a  predictor  which 

could  tolerate  impulse  responses  which  merely  became  constant  (rather 

than  zero)  for  large  T.  By  recomputing  ypzi-(0,t)  at  each  time  (n+1) 

for  which  a  prediction  is  needed,  an  expression  for  y  .j  can  be  obtained 

CL 

which  converges  whenever 

Hi_n  “  ht  =  0  (7.3) 

rather  than  the  usual  requirement  that 

H .  =  0  Vi  >T  (7.4) 

Specifically,  we  choose 

y»zi„ew(ntI,t)  ■X(H*"-H')u(t‘1,+ya(t) 


which  results  from  changing  the  upper  limit  in  equation  (7.2). 


7.3  Comparison  of  Predictors 

We  can  compare  the  behavior  of  these  two  predictors  by  considering 
the  actual  output  y2-(u,t)  resulting  from  zero  inputs  after  time  t.  We 
define  the  prediction  error  for  the  new  predictor  as 

e(n,t)  =  ypz1-(n,t)  -  y21-(n,t) 


Then 


T-n-1 

e(n+l,t)  =  J  (Hi+n  -  H.)u(t-i)  -  [yzi(n+l,t)  -  yfl(t)] 

We  note  that  if  H  from  the  internal  model  is  exactly  the  real  H,  then 
e(s,t)  =  0  for  all  s  as  long  as  H  satisfies  equation  (7.3).  Thus,  even  if 
H  does  not  become  zero  for  large  T,  good  prediction  is  possible.  By 
comparison,  the  error  in  the  regular  computation  would  be 

T-n-1  T-l 

e(n+l,t)  =  l  (H.  -H.)u(t-i)  -  [y  .(n+l,t) -y  (t)]  -  7  H.u(t-i) 

i=0  1  n  1  Z1  3  T-k-1  1 

T-l 

=  -  l  H.u(t-i) 

T-n-1  1 

which,  in  general,  is  not  zero. 

Similar  expressions  can  be  obtained  for  cases  of  bias  and  scale 
factor  errors  in  the  impulse  response.  If  the  internal  model  (H1)  is 
equal  to  the  real  impulse  response  plus  a  bias,  i.e.. 
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vr 


then  the  above  error  becomes  0  for  the  new  technique  and 


T-l 

e(n+l,t)  =  l  (H.  +  b)u(t-i ) 
T-n-1  1 


for  the  regular  predictor.  For  scale  factor  mismatch. 


h:  * 


qHi 


then 


T-n-1 

e(n+l,t)  =  (q-1)  l  (H.  -H.)u(t-i) 

i=0  1  n  1 

=  (q-D(yzi(n,t)  -ya(t)) 

for  the  new  technique,  which  results  in  an  error  which  grows  with  pre¬ 
diction  distance  (i.e.,  as  n  grows).  The  regular  predictor  has  this 
error  plus  another  term: 

T-l 

e(n+l,t)  =  (q-l)(y  ,(n,t) -y  (t))  -  £  qH.u(t-i) 

21  a  T-n-1  1 

In  general,  it  appears  that  the  modified  technique  is  better  under 
most  circumstances,  with  only  a  slightly  greater  computation  cost  (re¬ 
computing  many  "current  estimates"  based  on  prediction  lag).  This 
method  is  now  used  on  most  of  our  programs.  We  note  that  the  predictor 
for  n  steps  ahead  is  always  written  as  a  function  of  the  current  obser- 

ef 

vation,  and  not  as  a  function  of  the  n-1-  prediction.  This  preserves 
the  separation  between  prediction  and  control  computations,  since  the 
future  inputs  would  be  required  for  the  full  estimate  of  the  n-1^  term. 
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7.4  ARMA  Models 

The  use  of  past  inputs  and  impulse  responses  for  output  prediction 
corresponds  to  a  "moving  average"  filter  (see  Box  and  Jenkins  (1976)  for 
a  full  discussion  of  this  topic).  The  inclusion  of  the  current  observa¬ 
tion  creates  a  mixed,  autoregressive-moving  average  (ARMA)  predictor  as 
often  used  in  time  series  analysis.  It  was  noted  in  the  last  section 
that  the  inclusion  of  the  current  observation  permitted  the  treatment 
of  infinite-impulse  response  models.  Indeed,  this  is  a  familiar  fact  in 
ARMA  model ing,  and  often  motivates  the  use  of  mixed  models. 

In  our  current  version  of  IDCOM,  a  SO^1  order  impulse  response 
(moving  average)  is  used  with  a  first  order  autoregressive  term.  It 
seems  probable  that  a  lower  order  mixed-model  would  give  equally  good  (or 
better)  prediction  with  less  chance  of  numerical  problems,  straightforward 
identification,  and  the  ability  to  handle  measurement  noise.  Such  low- 
order  ARMA  models  combined  with  a  short  impulse  response  control  calcu¬ 
lation  would,  we  believe,  provide  good  control  and  estimation  for  systems 
in  noisy  environments  while  preserving  the  robustness  properties  of  MAC. 

We  intend  to  investigate  this  area  in  the  near  future. 

7.5  Optimal  Estimation 

Before  ending  this  section,  we  wish  to  briefly  outline  the  solution 
of  an  optimal  estimation  problem  for  IDCOM,  with  connections  to  ARMA 
models  and  simplified  estimation.  For  convenience,  we  begin  with  a  con¬ 
tinuous-time  plant  with  white  Gaussian  process  (w)  and  measurement  (v) 
noise,  and  dynamics: 
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x(t)  =  Ax  +  Bu  +  w(t) 
y  =  Cx 

z(t)  =  y(t)  +  v(t) 

where  x  is  the  state,  y  the  output,  and  u  a  known  input.  The  standard 
Kalman  filter  for  this  system  would  be 

X  =  AX  +  Bu  +  K(z-  Cx) 

«  (A-KC)x  +  Bu  +  Kz 

where  K  is  the  gain  matrix.  In  steady  state,  and  with  x(0) =  0,  this  system 
can  be  represented  as 


A 

X 


•t 

4>n  (t,i)[Bu  +  Kz]  di 
0  LL 


and 


9  = 


C4)^^(t,x)[Bu  +  Kz]  di 


where  <J>  is  the  closed-loop  transition  matrix 


♦CL(t,t)  •  e(A-KC){t-t) 


The  estimate  y  of  the  output  can  be  approximated  in  discrete-time  by 


N 

y(NA)  =  l  C<p~.  [Bu  +  Kz]A 
0  LL 
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which  can  be  written 


N  N 

y(NA)  =  l  Ht(i)uN  .  +  l  H?(i)zN  • 

1=0  1  N  i=0  £  N  1 

where  and  h^-C^^KA.  Without  noise,  K=0,  and  y  is  given  by 

y(t)  =  |  C4>0L(t,T)Bu(x)dT  +  y(0) 

where 

4>0L(t,T)  =  eA(t'x) 

and 

*  jQ  H0L(i)uN-i 

where  hqL  =  ^0LBA’  Thus,  we  note  the  similarity  between  noise-free, 

"open  loop"  prediction  based  on  u  alone  and  closed-loop  prediction  based 
on  both  inputs  and  measurements.  In  general,  there  is  no  "easy"  rela¬ 
tionship  between  Hq^  and  and  H^,  although  straightforward  calculations 
yield  the  latter  two. 

In  applying  these  results  to  IDCOM,  one  is  forced  to  create  50^ 
order  state  systems  to  handle  an  arbitrary  50^  order  impulse  response. 
The  calculations  for  such  large  systems  are  awkward  at  best,  and  the 
identification  of  lower-order  ARMA  models  would  greatly  simplify  the 
filter  computations. 
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Finally,  we  suspect  that  the  combination  of  exponentially-aged  fil¬ 
ters  with  reduced-order  ARMA  plant  models  may  prove  to  be  quite  workable 
and  practical.  Exponential  aging  (see  Miller,  1971,  or  Schweppe,  1973), 
sometimes  called  exponential  least  squares,  is  a  method  of  filter  design 
which  exponentially  weights  the  past  data.  By  introducing  an  exponential 
aging  term  on  the  measurement  noise  covariance,  the  Kalman  filter  gains 
remain  open  (non-zero)  and  the  filter  accepts  new  information  even  without 
a  process  noise  term. 

Currently,  the  construction  of  noise  covariances  is  often  fictitious, 
done  merely  to  keep  the  filter  gains  in  c.  region  which  results  in  good 
performance.  An  alternate  approach,  somewhat  less  arbitrary,  would  re¬ 
quire  control  system  designers  to  specify  only  measurement  noise  covar¬ 
iances  (measured  or  from  sensor  specifications)  and  an  exponential  weighting 
factor  (based  on  desired  bandwidth  but  explainable  as  easily  as  the 
reference  trajectory  time  constant  of  IDCOM).  Such  specifications,  along 
with  a  low  order  ARMA  process  model,  may  result  in  a  simple  estimator 
with  adequate  performance  and  the  ability  to  be  easily  modified  by  test 
personnel  as  sensors  change  or  experience  with  the  estimator  increases. 


7.6  Summary 

This  section  has  described  our  preliminary  work  in  improving  output 
prediction  for  IDCOM.  An  improved  predictor  was  described  for  normal 
operation  with  model  uncertainty.  Techniques  for  designing  predictors  in 
noisy  environments  were  discussed,  and  directions  for  future  research  were 


SECTION  VIII 


MAC  APPLICATION  TO  A  MISSILE  CONTROL  PROBLEM 


Overview 

This  section  describes  the  application  of  MAC  to  a  missile  attitude 
control  simulation.  The  control  program  for  MAC,  called  IDCOM,  is  dis¬ 
cussed  in  Section  6.1.  Two  control  computation  algorithms  used  in  the 
simulations  are  described  in  Sections  8.2  and  8.3.  The  basic  missile 
model  is  presented  in  Section  8.4,  and  modifications  made  to  the  roll 
axis  are  given  in  Section  8.5,  along  with  step  responses  for  the  modified 
(compensated)  missile.  All  simulation  results  will  be  presented  in 
Section  IX. 

8.1  IDCOM  Description 

MAC  is  generally  implemented  by  a  computer  program  called  IDCOM 
(for  IDentification  and  COftnand)  developed  by  ADERSA/GERBIOS  in  France. 

In  order  to  simplify  program  modifications  for  our  research,  another 
version  of  IDCOM  was  written  by  us  for  simulating  non-adaptive  MAC  appli¬ 
cations.  Our  initial  IDCOM  is  basically  equivalent  to  the  French  ver¬ 
sion,  with  some  modifications  as  described  below. 

A  block  diagram  of  the  components  of  our  IDCOM  is  shown  in  Figure  24. 
As  indicated  in  the  figure,  when  a  new  measurement  is  made,  it  is  fed 
to  two  blocks  of  IDCOM  and  used  to  compute  a  reference  trajectory  and 
a  "zero  input"  prediction  of  the  future  outputs  over  a  short  horizon  (for 
optimization).  The  reference  trajectory  is  a  first-order  exponential 
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Figure  24.  IDCOM  components. 


drawn  from  the  current  (measured)  output  to  a  given  set  point.  The  de¬ 
signer  supplies  a  time  constant  (x. )  for  this  exponential  for  each  output 
i.  The  "zero-input"  prediction  uses  the  past  inputs,  measurements,  and 
internal  impulse-response  model  to  predict  the  future  outputs  in  the  absence 
of  future  control  inputs. 

These  two  trajectories  (reference  and  zero-input  prediction)  are 
differenced  to  obtain  an  error  trajectory  to  be  minimized  by  the  future 
controls.  The  control  calculation  block  then  performs  this  minimization 
in  one  of  several  ways.  Once  an  input  sequence  has  been  computed,  the 


y,(t+l)> 
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first  input  is  applied  to  the  plant,  and  the  cycle  starts  again  after  the 
next  measurement. 

For  our  simulations,  the  zero-input  prediction  was  obtained  by  the 
"modified"  (bias-compensated)  estimator  discussed  in  Section  VII.  The 
controls  were  calculated  in  one  of  two  ways,  as  discussed  below. 

8.2  Basic  Control  Computation 

Our  first  control  computation  block  is  essentially  a  simple  inversion 
routine  (in  the  absence  of  constraints)  which  finds  inputs  to  zero  the 
error  between  the  reference  and  zero-input  trajectories.  To  make  the 
problem  tractable,  the  error  is  only  considered  at  a  few  points  in  the 
future,  and  the  future  controls  are  required  to  be  constant  over  inter¬ 
vals  ("blocks")  between  the  output-matching  points.  For  example,  three 
input  values  Up  u2  and  u3  may  be  computed  to  define  the  next  five  future 
controls  as: 

u(l,t)  =  u3 
u(2,t)  =  u2 
u(3,t)  =  u2 
u(4,t)  =  u3 
u(5,t)  =  u3 


where  u(s,t)  is  the  computed  control  at  time  t  to  be  applied  s  steps  in 
the  future.  This  is  shown  pictorially  in  Figure  25. 


current  12345 
time 


steps 

into 

future 


Figure  25.  Input  blocking. 

For  this  case,  the  outputs  to  be  controlled  (the  error  to  be  reduced) 
occur  at  the  endpoints  of  each  block,  i.e.,  1,  3  and  5  steps  in  the  fu¬ 
ture.  The  advantage  of  this  blocking  is  that  it  permits  reasonably  long 
optimization  horizons  with  low-dimensional  required  calculations.  In 
this  example,  only  three  numbers  are  computed  for  five  steps  ahead.  Since 
the  controls  will  generally  be  recalculated  one  step  ahead,  this  technique 
sacrifices  little. 

The  optimization  routine  used  in  the  basic  algorithm  tries  to  invert 
the  system  to  find  inputs  which  result  in  perfect  output  matching  at  the 
chosen  points.  If  this  is  possible,  then  the  first  control  causes  the 
first  output  to  be  correct  independently  of  the  future  inputs  and  outputs. 
If  constraints  are  encountered,  or  if  the  computation  time  is  shortened. 


8.3  Gradient  Algorithm 

In  order  to  improve  the  control  calculation  for  general  systems, 
a  gradient-projection  algorithm  was  substituted  for  the  basic  (inversion) 
control  computation.  This  new  algorithm  retains  the  input  blocking  dis¬ 
cussed  above  (to  reduce  the  dimension  of  the  problem)  but  does  not 
block  the  outputs.  Instead,  the  new  algorithm  minimizes  the  sum  of  the 
squares  of  the  output  errors  at  each  step  in  the  future  up  to  the  end 
of  the  last  input  block.  Thus,  for  the  example  above,  the  new  algorithm 
would  compute  three  control  numbers  to  minimize  five  future  errors. 

The  capability  also  exists  in  the  new  algorithm  for  adding  input 
and  output  weighting  matrices  if  desired.  In  general,  we  use  only  output 
weighting  unless  control  cost  is  meaningful.  Output  weighting  alone  is 
very  convenient,  permitting  easy  tuning  of  the  controller  to  achieve 
desired  performance  for  each  output.  Finally,  we  note  that,  for  only 
one  block  with  endpoint  one  step  ahead,  the  gradient  algorithm  and  in¬ 
version  routine  are  equivalent  (without  control  weights). 

8.4  Simulation  Model 

In  order  to  verify  the  theoretical  results  obtained  earlier  and  to 
demonstrate  the  behavior  of  MAC  in  an  aerospace  environment,  a  missile 
attitude  control  simulation  was  developed.  A  simple,  three-axis  attitude 
control  model  with  independent  pitch  axis  and  coupled  roll -yaw  dynamics 
was  chosen  from  AFIT  (1978).  The  model  represents  a  hypothetical  air- 


to-air  missile  with  asymmetric  aerodynamic  properties. 

The  model  has  six  states,  three  inputs  and  three  outputs,  with 


The  states  are: 


Xj  =  a  angle  of  attack  (rad) 

x2  =  q  perturbed  pitch  rate  (rad/s) 

x3  =  8  sideslip  angle  (rad) 

x^  =  p  perturbed  roll  rate  (rad/ s) 

x5  =  r  perturbed  yaw  rate  (rad/s) 

Xg  =  <p  roll  angle  (rad) 

with  inputs 

u1  =  elevator  angle  (6^)  (rad) 
u^  =  aileron  angle  (6p)  (rad) 
u3  =  rudder  angle  (6r)  (rad) 

outputs 

Yj  =  angle  of  attack  (a)  (rad) 
y2  =  sideslip  angle  (s)  (rad) 
y3  =  roll  angle  (<p)  (rad) 

and  parameters 

v  =  forward  speed 

Z,  =  dimensional  variation  of  z-force  with  downward  velocity,  sec 
w 

M  =  dimensional  variation  of  pitching  moment  with  angle  of 
a  attack,  sec"2 

Yft  =  dimensional  variation  of  y-force  with  sideslip  angle,  sec"* 
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Lfi  =  dimensional ^variation  of  rolling  moment  with  sideslip 
^  angle,  sec’^ 

Nr  =  dimensional?variation  of  yawing  moment  with  sideslip 
y  angle,  sec-<i 

M.  =  dimensional  variation  of^pitching  moment  with  pitch  control 
q  surface  deflection,  sec"^ 

Lo  =  dimensional  variation  of?rolling  moment  with  roll  control 
p  surface  deflection,  sec"^ 

Nxr  =  dimensional  variation  of^yawing  moment  with  yaw  control 
surface  deflection,  sec"^ 

g  =  acceleration  of  gravity 

<(>j  =  equilibrium  roll  angle 

Two  flight  conditions  were  chosen  for  study.  The  first  has  the  missile 
at  Mach  2  at  20,000  ft,  and  weighing  239.5  lb.  The  pitch  angle  (equi¬ 
librium)  is  9°  and  sideslip  is  0°.  The  second  flight  condition  is  Mach  4 
at  70,000  ft,  same  weight,  and  equilibrium  pitch  and  sideslip  of  14°  and 
11°  respectively. 

The  parameters  for  flight  condition  1  were: 

1  =  -1.4868 
w 

M  =  -149.93 

a 

M,  =  -281.11 
6q 

Yg  =  -.91237 
Lq  =  -1559.2 

P 

L.  =  8770.6 

6p 

Nq  =  290.48 

P 

N.  =  281.11 
or 
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The  parameters  for  flight  condition  2  were: 


Zw  =  -.27877 

M  =  -64.928 
a 

M(5q  =  -43.285 
Yg  =  -.25089 
Lg  -  1650.5 
Lr  =  1875.7 

6p 

N.  =  50.499 

p 

N6r  =  43.285 

The  models  are  considered  acceptable  for  variations  of  3°  in  a  and  0. 


8.5  Model  Modifications 

An  early  impulse  response  analysis  of  these  dynamics  indicated  a 
very  severe  roll  instability.  Since  IDCOM  works  best  with  finite  impulse 
responses,  we  chose  to  add  roll  angle  and  rate  feedback  to  the  aileron 
command,  thus  creating  a  compensated  system  for  IDCOM  to  control.  This 
changed  the  fourth  row  of  the  A  matrix  to  be: 


[0 


“LSpG<J> 


"L<5pV 


Due  to  the  roll -yaw  decoupling  from  the  pitch  axis,  this  change  did  not 
affect  the  angle-of-attack  dynamics  at  all. 

In  the  above  row,  Lgp  varies  with  flight  condition,  but  G^,  the 
compensator  gain,  was  fixed  to  be  h  in  our  tests.  AFIT  (1978)  presented 
values  of  of  0.127  at  flight  condition  1  and  1.056  at  flight  condition  2, 
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but  we  elected  to  use  a  simple,  fixed  gain  compensator  to  merely  stabilize 
the  plant.  IDCOM  would  be  used  to  get  the  required  performance,  and 
thus  an  "optimal"  should  not  be  needed. 

Step  responses  of  this  compensated  plant  to  control  inputs  are  shown 
in  Figures  26,  27  and  28  for  flight  condition  1  and  Figures  29,  30  and  31 
for  flight  condition  2.  A  sample  rate  of  10  Hz  (10  samples  per  second) 
was  used  for  these  plots,  with  linear  interpolation  (by  the  plotting  rou¬ 
tine)  filling  in  between  data  points.  As  shown  in  the  plots,  the  pitch 
axis  dynamics  are  quite  oscillatory  for  both  flight  conditions.  An 
analysis  of  the  decoupled  pitch  dynamics  revealed  that  flight  condition  1 
has  a  natural  frequency  of  12.24  r/s  (1.95  Hz)  and  a  damping  ratio  (c)  of 
0.061,  while  condition  2  has  a  frequency  of  8.06  r/s  (1.28  Hz)  and  a 
damping  ratio  of  0.017. 

8.6  Summary 

This  section  has  described  the  application  of  MAC  to  a  missile  atti¬ 
tude  control  simulation.  Two  versions  of  the  control  algorithm  (IDCOM) 
were  described,  and  their  differences  discussed.  The  missile  model  was 
presented,  and  changes  for  our  studies  were  noted.  The  next  section  will 
present  the  results  of  these  simulations. 
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Figure  26.  Step  response  to  elevator  input,  flight  condition  1. 
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Figure  27.  Step  response  to  aileron  input,  flight  condition  1. 
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Figure  29.  Step  response  to  elevator  input,  flight  condition  2. 
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Figure  30.  Step  response  to  aileron  input,  flight  condition  2. 
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A  =  angle  of  attack 
B  =  sideslip  angle 
C  =  roll  angle 
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Figure  31.  Step  response  to  rudder  input,  flight  condition  2. 
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SECTION  IX 


SIMULATION  RESULTS 


Overview 

This  section  presents  the  results  of  MAC  applied  to  a  missile  attitude 
control  simulation  as  described  in  Section  VIII.  The  first  section  (9.1) 
lists  the  simulation  parameters  and  baseline  conditions  for  the  majority 
of  the  runs.  Section  9.2  then  shows  the  performance  of  both  the  basic 
(inversion)  algorithm  and  the  modified  (gradient)  algorithm  for  the  base¬ 
line  conditions  and  some  slightly  different  parameter  values.  Changes  in 
the  reference  trajectory  time  constant  are  shown  in  Section  9.3,  and  the 
effects  of  input  rate  constraints  demonstrated  in  Section  9.4.  Section 
9.5  shows  the  behavior  of  the  algorithms  when  measurement  noise  is  pre¬ 
sent,  and  Section  9.6  demonstrates  robustness  limits  for  cases  where  the 
internal  model  and  real  missile  are  not  the  same. 

9.1  Simulation  Parameters 

The  following  sections  will  present  a  number  of  plots  showing  con¬ 
trol  responses  as  several  different  parameters  and  conditions  are  varied. 

In  order  to  facilitate  comparison  between  related  plots,  the  scales  have 
been  kept  constant,  if  possible,  within  each  series  of  runs. 

Unless  otherwise  noted,  the  following  conditions  existed  in  the  simu¬ 
lations: 

•The  sample  time  (A)  was  0.1  second. 
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•The  controls  were  computed  for  three  blocks  ending  at  one,  three, 
and  five  steps  in  the  future. 

•The  reference  trajectory  time  constant  (t)  was  0.1  seconds  (all 
axes) . 

•No  input  constraints  were  imposed. 

•The  set  points  were  changed  at  0.4  seconds  from  0°  to  15°  for 
angle  of  attack  (a)  and  10°  for  sideslip.  The  roll  set  point 
remained  at  0°. 

•For  the  gradient  algorithm,  the  output  weights  w^  were  all  equal 
(1);  no  input  weights  were  used. 

•Flight  condition  1  was  used. 

The  sampling  rate  of  10  Hz  corresponds  to  a  Nyquist  frequency  of 
5  Hz,  and  thus  all  plant  dynamics  should  be  slower  than  4  Hz.  We  expected 
the  closed  loop  plant  to  typically  be  controlled  with  a  0.1  second  ref¬ 
erence  trajectory  time  constant,  which  would  qualitatively  match  the 
performance  shown  in  AFIT  (1978).  This  time  constant  is  equal  to  a  1.59  Hz 
natural  frequency,  which  is  sufficiently  below  the  5  Hz  limit  for  these 
simple  simulations. 

For  the  simulation,  the  impulse  responses  for  IDCOM's  internal  model 
were  obtained  numerically  from  the  state  space  parameters.  The  simulation 
itself  ( i . e . ,  the  plant)  used  a  discrete-time  state  model  obtained  from 
solving  for  the  transition  matrix  for  the  state.  The  impulse  responses 
and  transition  matrix  calculations  were  carried  out  with  a  0.001  second 
step  size  for  accuracy.  This  method  insured  that  accurate  dynamic  be¬ 
havior  would  be  seen  at  the  sample  points,  although  no  attempt  was  made 
in  our  tests  to  examine  the  response  between  samples. 
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i  We  note  that  several  responses  are  outside  of  the  linear  range  for 

i 

'  the  model,  and  thus  should  not  be  considered  as  absolute  results.  Rather, 

j 

i  the  inputs  and  outputs  should  all  be  scaled  down  (to  less  than  3°)  if 

absolute  numbers  are  desired.  Since  all  of  the  operations  (except  for 

;  the  constraint  case)  are  linear,  this  presents  no  problems. 

-  We  also  note  that  the  low  pitch  damping  ratio  of  the  missile  at 

'  each  flight  condition  presents  a  difficult  control  task  for  any  digital 

controller.  An  actual  missile  implementation  would  benefit  from  an  analog 
stabilization  loop  inside  of  the  IDCOM  control  loop  (as  was  modeled  for 
the  roll  axis  and  discussed  in  Section  VIII).  We  did  not  use  this  internal 
control  structure  for  angle  of  attack,  however,  since  we  wanted  to  ex¬ 
plore  the  operation  and  limitations  of  different  IDCOM  techniques  in  this 
environment. 

*•  9.2  Initial  Control  Tests 

This  series  of  tests  demonstrates  the  basic  control  behavior  of  the 

,  14 

two  versions  of  IDCOM  tested.  The  first  figure  (Figure  32)  shows  the 
basic  (inversion)  algorithm  for  the  baseline  conditions  described  above. 
Figure  33  shows  the  same  controller  but  only  looking  one  step  ahead 
(i.e.,  one  endpoint,  one  step  ahead).  We  see  that  without  constraints, 
the  results  are  nearly  identical.  The  large  oscillations  in  inputs  are 
due  to  the  inversion  of  an  oscillatory  system  (the  missile).  The  con¬ 
trolled  outputs  are  very  good,  but  the  control  behavior  appears  unde¬ 
sirable. 

^All  figures  for  Section  IX  are  grouped  at  the  end  of  the  section, 
beginning  on  page  141. 
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The  next  five  plots  show  the  modified,  gradient  algorithm.  Figure  34 
shows  the  baseline  conditions,  and  Figure  35  shows  the  response  to  a  10° 
roll  command  with  0®  sideslip.  Note  the  much  smoother  control  sequence 
compared  to  the  inversion  routine  (Figure  32).  Figure  36  shows  the  gra¬ 
dient  algorithm  only  looking  one  step  ahead,  and  thus  converging  to  the 
inversion  technique.  Figure  37  shows  the  improvement  of  just  considering 
two  endpoints  (one  and  five  steps  ahead),  and  Figure  38  considers  three 
endpoints  (one,  three,  and  10  steps  ahead),  using  a  longer  optimization 
horizon.  This  results  in  slightly  slower  initial  response  but  smoother 
controls  than  Figure  34. 

9.3  Reference  Trajectory  Changes 

These  tests  show  the  effect  of  reference  trajectory  time  constant  on 
speed  of  response.  The  gradient  algorithm  was  used  for  this  series  and 
results  should  be  compared  to  Figure  34.  Figure  39  shows  the  response 
for  a  time  constant  of  .5  second  (from  .1  seconds  in  Figure  34),  while 
Figure  40  speeds  the  reference  to  0.05  seconds.  For  comparison,  the 
fast  time  constant  (.05)  with  a  one-step  look  ahead  controller  (gradient 
here,  but  equivalent  to  the  inverse)  is  shown  in  Figure  41. 

9.4  Input  Rate  Constraints 

This  series  of  runs  shows  the  effect  of  adding  rate  constraints  to 
the  control  inputs.  IDC0M  directly  considers  magnitude  and  rate  constraints 
in  its  control  calculation  so  that  improved  compensation  is  possible  (i.e., 
other  controls  may  be  adjusted)  when  one  input  is  limited.  For  the 
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missile  control  simulation,  magnitude  limits  did  not  seem  important, 
since  most  control  surfaces  had  large  travel  available.  Rate  limits 
(i.e.,  derivative  of  control  surface  motion)  did  seem  realistic,  however, 
and  we  attempted  to  investigate  their  effects  on  performance. 

We  recall  that  the  baseline  tests  showed  large  swings  in  control 
inputs  were  needed  for  the  inversion  algorithm  to  achieve  its  excellent 
tracking.  These  input  rates,  although  alarming,  barely  exceed  a  250°/ 
second  rate  limit  suggested  by  the  Air  Force  for  this  simulation.  Fig¬ 
ure  42  shows  the  response  of  the  inversion  algorithm  (one  step  look 
ahead)  to- a  rate  limit  of  250°/second.  The  performance  is  generally 
indistinguishable  from  that  of  the  baseline  test.  Figure  32.  A  very 
strict  50°/second  rate  limit  is  imposed  in  Figure  43,  and  this  time  the 
output  does  not  perfectly  match  the  reference  trajectory  but  the  control 
swings  are  nearly  eliminated. 

9.5  Effects  of  Noise 

This  series  of  tests  demonstrates  the  behavior  of  IDCOM  in  the  pre¬ 
sence  of  white,  Gaussian  measurement  noise.  The  gradient  controller  was 
used  for  these  runs,  although  the  first  example  (Figure  44)  uses  only 
one  optimization  point  and  is  the  equivalent  of  the  inversion  technique. 
The  controller  was  used  in  a  regulator  mode,  with  no  set  point  step  com¬ 
mands,  and  the  plot  scale  was  reduced  to  show  the  output  noise.  The 
plots  are  actually  of  the  measured  output  including  the  true  plant  out¬ 
put  plus  measurement  noise.  They  thus  show  slightly  more  motion  than 
the  true  plant  output  alone  exhibited.  The  noise  was  chosen  to  have  a 
variance  of  0.25,  applied  at  the  sample  time,  0.1  seconds  apart. 
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Figure  44  shows  the  outputs  for  the  inversion-equivalent  case.  The 
actual  output  noise  variances  were  0.1348,  0.0922,  and  0.1001  for  angle 
of  attack,  sideslip,  and  roll  respectively.  Figure  45  shows  the  gradient 
algorithm  looking  three  steps  ahead  (as  in  the  baseline  case).  The  output 
noise  variances  for  this  case  are  0.1016,  0.0387,  and  0.0290  for  angles 
of  attack,  sideslip  and  roll,  respectively. 

Although  it  is  dangerous  to  draw  conclusions  from  too  small  a  random 
sample,  the  gradient  algorithm  appears  to  retain  its  smooth  control 
properties  (compared  to  the  inversion  technique)  with  little  or  no  penalty 
in  output  noise  variance.  The  much  more  easily  analyzed  inversion  rou¬ 
tine  may,  therefore,  provide  a  rough  estimate  of  noise  behavior  for  the 
gradient  technique. 

9.6  Model  Mismatch 

This  section  demonstrates  the  behavior  of  the  IDC0M  algorithms  when 
the  internal  (to  IDC0M)  impulse  response  model  is  not  the  same  as  the 
state  space  plant  model  used  to  simulate  the  missile.  Two  different 
types  of  mismatch  are  shown  below:  a  simple  gain  error  as  discussed  in 
Section  3.1,  and  a  flight  condition  mismatch  (e.g.,  internal  model  of 
flight  condition  1  with  plant  from  flight  condition  2)  as  treated  in 
Appendix  B. 

9.6.1  Gain  Mismatch 

As  discussed  in  Section  3.1,  when  the  internal  model  impulse  response 
is  equal  to  a  constant  gain  (1/q)  times  the  real  impulse  response,  the 
controlled  system  is  stable  if 


(9.1) 


0 


where 


-A/t 

a  =  e 

for  sample  time  A  and  reference  trajectory  time  constant  r.  This  assumes 
that  a  stable  system  is  being  controlled,  and  that  the  inversion  algorithm 
is  being  used. 

To  demonstrate  this  result  we  used  the  gradient  controller  with  one 
endpoint  one  step  ahead  (equivalent  to  the  inversion  routine)  and  set 
point  step  commands  of  3°  (angle  of  attack),.  2°  (sideslip)  and  0°  (roll). 
These  set  points  have  the  same  relative  value  as  the  baseline  conditions 
an«j  permitted  use  of  the  same  plot  scaling. 

For  the  normal  sample  time  (A  =  0.1  seconds)  and  time  constant  (t  =  0.1 
seconds),  equation  (9.1)  indicates  instability  for 

q  >  3.164 

Figure  46  shows  the  case  where  q  =  3,  and  indeed  the  system  is  barely  stable. 
Increasing  q  to  4  resulted  in  instability  (not  shown). 

As  indicated  by  equation  (9.1)  for  any  given  q,  robustness  can  be 
improved  by  increasing  the  reference  trajectory  time  constant.  For  q  =  4, 

r  >  0.144  seconds 

should  produce  stability.  Indeed,  a  case  with  t  =  0.15  is  shown  in  Fig¬ 
ure  47,  and  the  system  is  again  stable,  but  not  by  much.  Increasing 
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r  to  0.2,  as  shown  in  Figure  48  (q  =  4),  made  the  system  much  more  stable 
(less  oscillatory) . 

Although  no  theoretical  results  exist  yet  for  the  gradient  routine, 
its  generally  smoother  control  was  expected  to  result  in  greater  robust¬ 
ness  than  the  inversion  algorithm.  This  proved  to  be  true  in  our  limited 
testing.  Figure  49  shows  the  gradient  algorithm  with  three  endpoints 
(baseline  conditions),  q  =  3,  and  t  =  0.1  seconds  as  in  Figure  46.  The 
gradient  algorithm's  behavior  is  clearly  superior.  Indeed,  for  q  =  4, 
t  =  0.1  seconds,  when  the  inversion  case  was  unstable,  the  gradient  al¬ 
gorithm  gave  reasonably  good  performance,  as  shown  in  Figure  50. 

These  tests  seem  to  both  confirm  the  theory  of  Section  3.1  and  imply 
that  the  robustness  bounds  for  the  simple  inversion  are  conservative 
limits  for  more  sophisticated  optimization  algorithms. 

9.6.2  Flight  Condition  Mismatch 

One  of  the  most  attractive  possibilities  of  robust  control  of  air¬ 
craft  and  missiles  is  the  potential  for  reducing  the  amount  of  gain 
scheduling  required  as  flight  conditions  change.  Robustness  analysis 
for  this  case  can  be  quite  tedious,  as  indicated  in  Appendix  B,  and  one 
often  must  simulate  and  examine  controllers  to  determine  the  best  struc¬ 
ture  and  gains.  This  section  presents  a  very  limited  demonstration  of 
the  potential  of  IDCOM  in  these  areas. 

The  two  flight  conditions  chosen  (see  Section  8.4)  represent  widely 
different  systems,  each  difficult  to  control.  Figures  34  and  36  earlier 
showed  the  performance  of  IDCOM  (gradient  version)  for  three  and  one 
endpoint  optimization  (respectively)  at  flight  condition  1.  Similar 


138 


baseline  plots  for  flight  condition  2  are  given  in  Figures  51  (three 
endpoints)  and  52  (one  endpoint),  but  with  step  commands  of  3°  (angle 
of  attack),  2°  (sideslip)  and  0°  (roll).  Figure  52  shows  the  near  in¬ 
stability  of  the  inverse  control  technique  for  the  highly  undamped 
pitch  dynamics  at  this  flight  condition.  Again,  the  gradient  technique 
(Figure  51)  is  much  better.  In  all  of  these  runs,  the  internal  model 
was  correctly  identified  (off-line)  for  each  flight  condition. 

We  next  tried  mixing  the  two  internal  models  and  actual  flight  con¬ 
ditions,  to  see  if  one  control  model  would  work  at  the  two  extreme 
conditions.  We  simulated  both  gradient  and  inversion  techniques,  and 
neither  was  stable  with  the  mismatched  model  for  a  reference  x  of  0.1 
seconds.  We  then  tried  slowing  the  reference  trajectory  down,  and  were 
able  to  stabilize  the  systems,  with  the  consequent  sluggish  performance. 

Figure  53  shows  the  inversion  (actually  gradient  with  one  endpoint) 
algorithm  (x  =  l  second)  at  flight  condition  2  and  an  internal  model  from 
flight  condition  1.  The  set  points  were  the  baseline  15°,  10°,  and  0°. 
Figure  54  shows  the  same  case  for  a  three  endpoint,  gradient  algorithm. 
The  responses  are  very  similar,  with  less  control  oscillation  in  the 
gradient  case. 

Figure  55  shows  the  inversion  algorithm  for  the  reverse  case  (model 
from  condition  2  and  missile  simulated  at  condition  1)  with  lower  set 
points  (3°,  2°,  and  0°).  The  missile  immediately  approaches  the  set 
point  because  of  the  unexpectedly  large  response  (compared  to  the  control 
calculation  model),  even  though  the  time  constant  is  still  1  second 
(reduced  for  stability).  Figure  56  shows  the  gradient  (three  point)  al¬ 
gorithm  for  this  case  with  the  familiar  improvement  in  control  smoothness 


9.7  Summary 


This  section  has  presented  the  simulation  results  from  applying  MAC 
to  a  missile  attitude  control  model.  Some  of  the  theoretical  results  of 
the  earlier  analysis  sections  have  been  demonstrated  for  a  simple  al¬ 
gorithm  which  approaches  the  idealized  MAC  considered  in  Section  III. 
Results  were  also  presented  for  an  improved  controller  currently  under 
development. 
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Figure  34a.  Gradient  algorithm  at  baseline  conditions. 
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Figure  34b.  Gradient  algorithm  at  baseline  conditions 
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Figure  36b.  Gradient  algorithm  with  one  endpoint  (becomes 
inversion  algorithm,  cf.  Figure  33). 
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Figure  40a.  Gradient  algorithm,  fast  trajectory. 
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Figure  40b.  Gradient  algorithm,  fast  trajectory 
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Figure  42a.  Inversion  algorithm,  250°/second  rate  limit 
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Figure  44a.  Inversion  algorithm  with  measurement  noise. 
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Figure  45b.  Gradient  algorithm  with  measurement  noise. 
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Figure  46a.  Inversion  algorithm  with  gain  mismatch  of  3,  t=  .1. 
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Figure  47a.  Inversion  algorithm  with  gain  mismatch  of  4,  T 
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Figure  50a.  Gradient  algorithm  with  gain  mismatch  of  4,  x=  .1. 
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Figure  51a.  Gradient  algorithm,  flight  condition  2  baseline. 


179 


INPUT  (degrees) 

-20*00  -10.00  0.00  10.00 


I - , - 1 - 1 - 1 - 1 — 

0.00  1-00  2.00  3.00  4.00  5-00 

TIME  (seconds) 


Aileron  Deflection 


Rudder  Deflection 


I - 1 - 1 - 1 - 1 - 1 — 

0.00  1.00  2.00  3.00  4.00  5-00 

TIME  (seconds) 


Figure  51b.  Gradient  algorithm,  flight  condition  2  baseline. 
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Figure  53a.  Inversion  algorithm,  control  model  from  condition  1 
with  missile  at  condition  2. 
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Figure  53b.  Inversion  algorithm,  control  model  from  condition  1 
with  missile  at  condition  2. 
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Figure  54a.  Gradient  algorithm,  control  model  from  condition  1 
with  missile  at  condition  2. 
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Figure  56a.  Gradient  algorithm,  control  model  from  condition  2 
with  missile  at  condition  1. 
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SECTION  X 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  overall  conclusion  of  this  study  is  that  MAC  technique  has  a 
sound  mathematical  and  empirical  basis.  It  is  a  highly  flexible,  in¬ 
tuitive,  and  general  approach  to  control  design  which  fully  exploits 
the  capabilities  of  current  microprocessors.  MAC  can  handle  hard  con¬ 
straints,  time  varying  system  characteristics,  and  unequal  numbers  of 
inputs  and  outputs.  These  features,  ease  of  implementation  and  its 
theoretical  properties  of  robustness  make  it  a  very  powerful  technique 
for  guidance  and  control  of  flight  vehicles. 

Specific  conclusions  of  this  study  are:  (i)  the  mathematical  pro¬ 
perties  of  an  idealized  MAC  can  be  analyzed  by  conventional  control  analysis 
techniques.  For  unconstrained  single-input,  single-output  minimum  phase 
systems,  MAC  design  is  similar  to  inverse  control.  An  appropriate  choice 
of  reference  trajectory  is  required  to  achieve  desired  robustness, 
tracking,  and  disturbance  rejection  properties. 

(ii)  Simple  analytical  criteria  can  be  derived  for  the  robustness  of 
MAC  in  terms  of  gain  margin,  which  agree  quite  well  with  the  simulation 
results.  It  is  shown  that  there  is  a  direct  relationship  between  MAC 
robustness  and  speed  of  response  of  the  reference  trajectory  (or  of  the 
closed  loop  system). 

(iii)  In  general,  the  performance  of  the  closed  loop  prediction  MAC 
is  better  than  that  of  the  open  loop  prediction  MAC.  However,  there  are 
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cases  (e.g.,  large  initial  distance  from  the  set  point)  where  a  combined 
use  of  open  loop  prediction  and  closed  loop  prediction  strategies  may  be 
useful . 

(iv)  For  nonminimum  phase  systems,  inverse  control  leads  to  unbounded 
inputs.  The  MAC  design  approach  can  be  easily  modified  to  overcome  this 
difficulty,  however.  The  modified  MAC  design  uses  two  different  models: 

a  prediction  model  which  is  a  close  replica  of  the  true  system  and  a  con¬ 
trol  (or  deconvolution)  model  which  has  a  stable  inverse.  An  optimal 
solution  to  the  latter  problem  is  given  in  terms  of  a  Riccati  equation. 
(The  use  of  the  Riccati  equation  solution,  in  this  context,  is  very  dif¬ 
ferent  from  that  in  LQR  design.)  The  optimal  design  performs  much  better 
than  pole  placement  and  weighting  techniques. 

(v)  MAC  design  can  be  extended  to  continuous  time  systems  and  the 
optimal  inputs  can  be  computed  by  parameter  optimization  methods  using 
suitable  basis  functions.  The  control  law  is  shown  to  be  a  piecewise 
linear  function  of  the  reference  trajectory  and  current  state  (or  equi¬ 
valently  the  past  inputs  and  outputs)  of  the  system. 

(vi)  For  continuous  time  systems  with  discrete  observations,  un¬ 
observed  outputs,  process  noise,  and  measurement  noise,  the  MAC  approach 
can  be  applied  in  conjunction  with  a  state  estimator  (Kalman  filter). 

The  control  computations  are  still  performed  using  an  impulse  response 
model  and  a  quadratic  programming  algorithm.  The  sample  rate  selection 
problem,  however,  can  be  studied  by  an  extension  of  the  LQG  theory.  The 
effect  of  sample  rate  on  MAC  performance  can  be  expressed  in  terms  of 
the  solution  to  the  well-known  estimation  and  control  Riccati  equations. 
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( vi 1 )  The  missile  control  simulations  demonstrated  superior  perfor¬ 
mance  achievable  with  MAC  in  a  missile  control  environment.  The  basic 
MAC  was  shown  to  be  tolerant  of  model  mismatch  and  measurement  noise, 
as  indicated  in  the  theoretical  development.  An  alternate  optimization 
algorithm  was  shown  to  preserve  the  good  output  tracking  of  the  regular 
algorithm  while  using  smoother  control  sequences.  The  alternate  algorithm 
also  appears  to  be  more  robust  than  the  basic  one,  and  just  as  tolerant 
of  noise. 

(viii)  Control  design  involves  subtle  tradeoffs  between  several  con¬ 
flicting  objectives  such  as  tracking,  robustness,  constraint  satisfaction 
and  noise  performance.  Using  the  MAC  approach,  a  large  number  of  designs 
can  be  tried  rapidly  and  the  above  tradeoff  can  be  made  easily.  MAC  de¬ 
sign  is  also  self-suggestive  for  improvements  since  it  reveals  the  factor 
limiting  performance. 

The  following  reconmendations  are  made  for  further  research:  (i) 
adaptive  MAC:  the  MAC  approach  is  easily  extended  to  the  case  where  the 
plant  model  is  varying  in  an  a  priori  unpredictable  fashion.  It  has  been 
shown  that  identification  and  control  are  dual  problems  in  the  MAC  formu¬ 
lation.  However,  the  theoretical  properties  of  adaptive  MAC  have  not  been 
studied  thoroughly.  It  is  recommended  that  the  convergence  and  robust¬ 
ness  properties  of  adaptive  MAC  be  studied  analytically  and  the  role  of 
test  signals  (explicit  and  implicit)  be  quantified  in  terms  of  control 
performance. 

(ii)  Impulse  response  representation  is  not  suitable  for  lightly 
damped  and  unstable  systems.  A  number  of  heuristic  techniques,  including 


compensation  prior  to  control  by  MAC  are  available  and  have  been  used  in 
practical  applications.  It  is  recommended  that  these  techniques  be  analyzed 
theoretically  to  select  optimum  designs  for  specific  situations. 

(iii)  Since  the  MAC  approach  is  equally  applicable  to  guidance  design, 
it  is  recommended  that  both  guidance  and  control  problems  for  flight 
vehicles  be  investigated  using  the  MAC  approach.  In  particular,  the  ap¬ 
plication  of  MAC  to  cruise  missile  guidance  is  very  promising  (Reid  et  al . , 
1980). 

(iv)  Several  algorithmic  improvements  in  MAC  design  have  been  described 
in  this  report.  It  is  recommended  that  these  improvements  be  tested  on 
more  complex  simulations. 

(v)  The  next  logical  step  in  MAC  development  is  testing  a  research 
aircraft  such  as  TIFS  and  NAVION.  This  would  involve  determination  of 
real-time  computer  requirements,  comparison  with  other  control  techniques 
and  performance  evaluation. 
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APPENDIX  A 


GENERAL  STUDY  OF  SINGLE  INPUT  SINGLE  OUTPUT  LINEAR  TIME 
INVARIANT  CONTROL  LAWS.  APPLICATION  TO  AN 
ADAPTED  MODELS  ALGORITHM  CONTROL  (AMAC) 


L.  PRALY 


Abstract:  In  this  study,  we  come  back  on  some  characteristics  of  linear 
time  invariant  control  laws  and  we  show  how  the  single  input  single  output 
(SISO)  adapted  model  algorithm  control  (AMAC)  is  a  technique  for  designing 
such  a  law. 
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This  appendix  (pp.  197-250)  originally  appeared  as  an  internal  report  at 
ADERSA/GERBIOS,  and  it  is  reproduced  here  in  its  original  form  with  its 
own  references  and  appendices. 


A  -  RETURN  ON  LINEAR  TIME  INVARIANT  LAWS  OF  CONTROL 


A- 1  Introduction 


Let  us  return  on  what  is  the  problem  of  sampled-data 
control  systems  synthesis. 

Given  a  mathematical  model  of  a  process,  functionnal 
operator  between  an  input  and  an  output,  given  a  set  of 
specifications,  given  a  method  to  compute  a  control  input, 
the  problem  of  synthesis  may  be  defined  as  follows  :  find 
the  parameters  used  in  the  computation  of  the  control  input 
such  that  the  mathematical  process  with  that  input  meets 
all  the  specifications.  We  spoke  about  a  mathematical  model 
of  a  process  and  not  about  a  real  physical  process.  We  will 
say  a  law  of  command  to  be  robust  if  it  can  be  used  on  a 
physical  process. 

More  precisely  we  call  robustness  the  coherence  between 
approximations  of  a  mathematical  representation  of  a  physical 
process,  and  the  sensitivity  of  performance  criteria  defined 
by  the  specifications,  to  variations  of  this  representation. 

Let  Pp  be  the  nominal  mathematical  process,  the  control 
input  is  designed  for,  if  the  performance  criteria  are  continuous 
in  Pq,  we  can  expect  the  satisfaction  of  the  specifications 
for  any  P  in  the  vicinity  of  P^.  Then,  the  physical  process 
we  want  to  command  must  have  representations  each  in  this 
vicinity. 

Another  way  to  formulate  the  problem  is  :  the  set  of 
mathematical  processes,  images  of  the  physical  process,  must 
be  enclosed  in  the  set  of  mathematical  p roc esses  which  verify 
the  specifications  for  a  given  control  law. 
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So  we  will  successively 


-  define  a  mathematical  process 

-  define  a  set  of  specifications 

-  define  a  control  law 

-  find  relationships  between  parameters  of  the  control 
law 

-  study  the  sensitivity  of  the  performance  criteria. 
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A- 2  Definitions 


A-2.1  Definition  of  the  mathematical  model  of  a  process 

A-2.1.1  Def inition 

We  define  here  a  discrete  time  mathematical  model 
of  a  process  as  a  transformation  of  a  set  of  sequences  called 
inputs  into  another  set  of  sequences  called  ouputs . 

We  differentiate  three  types  of  signals  between 
the  input  sequences  : 


noted  e 

n 


a  controlable  measured  signal  called  control  and 


an  uncontrolable  but  measured  signal  called 
measured  disturbance,  and  noted  v 

n 

an  uncontrolable  unmeasured  signal  called 

disturbance  and  noted  w  . 

n 

For  single  input  single  output  systems  e,v,w 
are  scalars  and  so  is  the  output  noted  s^. 

So,  if  P  is  an  operator  on  sequences,  we  have 
the  relation  between  inputs  and  output  : 


s ( . ) =P (e { . ) , v ( . ) ,w( . ) ) . 


A-2.1  .2  Hypotheses 

A - 2 . 1 . 2 . 1  Hypothesis  on  P(Hl)  : 

We  suppose  P  to  be  a  linear  time 
invariant  operator  which  is  of  rational  type  and  asymptoti¬ 
cally  stable.  Moreover  we  suppose  a  non  zero  static  gain. 
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A - 2 . 1 . 2 . 2  Hypothesis  on  disturbance s  (H2) 

We  suppose  both  measured  and  unmeasured 
disturbances  to  be  causal  and  to  admit  z  -  transforms  which 
verify  conditions  of  final  value  therorem  CO- 

If  the  disturbances  are  represented  as 
stochastic  processes,  these  hypotheses  are  made  on  the 
mathematical  expectations  and  all  the  following  deterministic 
results  must  be  considered  in  mathematical  expectation. 


Moreover  we  suppose  the  output  to  be 
lineraly  time-invariant  dependant  on  the  disturbances.  So  we 
introduce  a  new  linear,  time  invariant,  asymptotically  stable 
operator  Q  between  the  measured  disturbance  and  the  output. 


A- 2. 1.3  Representation  of  the  mathematical  model 
of  the  process 

With  hypotheses  Hi,  H2,  we  compute  the 
output  s(n),  from  the  inputs  e ( n ) , v ( n ) , w ( n)  by  the  recursive 
equation  : 


N, 


E  h  9i  en-l*  E  ”1  V"-l 

i=0 


N 


i  =  0 


1  =  0 


+  V  f  ,  w 
6*  1  n_1 


(1  ) 


where  ( f ± > L t(0  .  N  ) '  (gi)ic(0,.,N  )' 

f  9 


(h.)  .  ,  are  time  invariant  scalars 

l  i€  (0 ,  .  ,Nh> 


Neglecting  the  initial  conditions 
(justified  by  asymptotic  stability),  we  can  represent  (1) 
in  a  more  concise  way  using  z  -  transforms 


s(z)=P(z)e(z)+Q(z)v(z)+w(z) 


(2) 


where  : 


P  { z) 


Pn(z) 

Pd(z) 


n 

PnU)=  £_  gi 


N  -i 

„g 


i=0 


qnU) 


clzl'i7^T  • 


.  h  V1 

qn(zl"  2_  "l  2  1 

i  =  0 


pd(2) ■qd(z) 


-t 

i  =  0 


Nf-i 


f .  Z 
1 


( 


and  from  the  causality  principle,  degree  of  p^lq^)  is  ^ess 
than  degree  of  P^^d^’ 

Moreover  from  the  hypotheses,  the  roots  of  p^(z) 
and  q^tz)  are  strictly  in  the  unit  circle. 

So  we  get  the  block  representation  given  by 

figure  1  : 


FIGURE  1 


Representation  of  the  process 


-2.2  Definition  of  a  set  of  specifications 


A -2. 2.1  Output  regulation 

We  want  the  effect  of  non  decreasing  distur¬ 
bances  on  the  process  to  be, in  some  sense,  minimized  or 
e 1 iminated . 

A- 2 . 2 . 2  Output  tracking 

Given  an  external  non  diminishing  signal 

called  the  set  point  and  noted  u^  ,  we  want  the  output 

s  to  track  u  with  minimal  or  ideally,  zero  steady 
n  n 

state  error.  For  this  problem,  we  impose  a  causal  set  point 
with  z  -  transform  which  verifies  conditions  of  the  final 
value  theorem. 


A -  2 . 2 . 3  Internal  stability 

In  both  cases  it  is  also  imperative  that 
an  appropriate  control  law  be  designed  in  such  a  way  as  to 
insure  an  asymptotically  stable  design  i.e.  the  relations 
between  the  external  signals  (set  point,  measured  and 
unmeasured  inputs)  and  the  internal  signals  (control , output) 
must  be  stable  in  some  sense. 

A -2. 2. 4  Asymptotic  convergence 

We  will  summarize  the  preceding  definitions 
by  the  asymptotic  convergence  of  the  output  s^  to  the 

set  point  u 

n 


limu-s=0  (4) 

n  -•  oc  n  n 
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In  fact  we  have  only  here  the  least  constraints 
to  set  any  system  of  control.  The  synthesis  of  such  a  system 
must  also  take  into  account  the  behaviour  of  this  convergence 
and  need  performance  criteria  [2] .  From  greater  variability  we 
keep  ourselves  within  the  convergence  criterion. 


Y 


A-2.3  Definition  of  the  control  law 


A- 2. 3.1  Definition 


We  call  control  law  a  method  to  compute  future  controls 

given  the  observation  of  all  the  measurable  past  signals.  To 

get  a  very  general  linear  time-invariant  control,  we  compute 

a  future  control  en+i  *  given  the  past  measured  signals 

(e  ,  s  ,  u  ,  v  ;  m  <  n)  as  a  finite  linear  combination  : 

m  m  m  m  ' 


N 

e  (n+1  )  =  -\ 

i  =  0 


1  =  0 


v 

n-i 


(5) 


4 


> 

1*' 


Or  using  z  -  transforms,  we  write  : 

c(z)e(z)=r(z)u(z)-d(z)s(z)-b(z)v(z)  (6) 

with  c(z) ,  r(z) ,  d(z) ,  b(z)  z  -  polynomials  such  that 
degree  of  c(z)  is  greater  than  degree  of  r(z),  d(z)  or  b(z) 
and  c(z)  is  mutually  prime  with  r(z) ,  d(z)  and  b(z) . 

Note  that  from  the  homogeneity  of  equation  (6) ,  there 
is  no  use  to  take  rational  functions  instead  of  polynomials. 

A -  2  .  3  .  2  I nterpretation 

Equation  ((>)  has  the  block  diagram  representation  given 
in  figure  2. 


sa 
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and  the  regulation  feedback  and  feedforward  transfers  : 


s  t _c  (z)  Q  (z)  -b  (z)  P  (z) 
rv  v  *  c ( z) +d ( z) P ( z) 


(11) 


srw<’> 


_ c  (z) 

%  (z)  +d  (z)  P  (z) 


(12) 


.  b(z)+d(z)Q(z) 

rv  '  ”c(z)+d(z)P(z) 


(13) 


Er  (z)=- 


d  ( z ) 


•  w 


c(z)+d(z)P(z) 


(14) 


We  can  see  that  the  poles  of  any  transfer  are 
given  by  the  roots  of  the  expression  c  (z) +d (z) P  (z)  .  Moreover 
from  the  stability  of  P(z),  Q(z)  and  the  hypothesis  of 
mutual  primeness,  a  necessary  and  sufficient  condition  of 
internal  stability  is  given  by  the  stability  of  the  control 
and  more  precisely  by  the  stability  of  the  Efl(z)  transfer. 


We  shall  note  that  given  the  stability  conditions, 

the  sensor  d(z)  determine  the  Er  (z)  regulation  transfer,  the 

1  w 

compensator  c(z)  determines  the  S r^lz)  regulation  transfer 
and  r(z)  determines  the  E  (z)  tracking  transfer.  With 

A 

the  error  tracking  transfer  : 


1 


S 

a 


(z)  =S 


•  w 


.  ,  ,  <d (z) -r (z) )P (z) 
1  '  c(z)+d(z)P(z) 


(15) 


we  remark  that  the  difference  between  d(z)  and  r(z) 
d i f f er enci a te s  between  regulation  and  tracking  behaviours. 


Whith  expression  (11),  if  it  is  possible  to  get  : 

c(z)Q(z)-b(z)P(z)  (16) 

we  will  be  able  to  compensate  completely  the  measured  distur¬ 
bance. 
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From  expressions  (9),  (11),  (12),  we  can  write 

expression  (4)  using  the  final  value  theorem  : 


lim  (l-z) (  s (i)-u(z))«0 
z  -*  1  . 


(17) 


Supposing  s(z)>u(z)  to  be  defined  in  the  ring 
(l,+00  ).  in  order  to  separate  set-point,  measured  and 
unmeasured  disturbances  actions,  we  will  transpose  the  set 
of  specifications  into  four  constraints  : 

regulation  constraints  : 


Sr  (1)« 
rw 


c  (  1  ) 


c(l)+d(l)P(l) 


=  0 


(18) 


s  (1)  ,c(l)Q(l)rb(l)P(l)  = 
Srvu)'  c  ( 1 )  +d  ( 1 )  P  ( 1 )  ° 


(19) 


tracking  constraint  : 


S  (1)=1= _ L.(J2IJUL  _ 

a'  '  c(l)+d(l)P(l) 


(20) 


And  stability  constraint 


The  roots  of  c(z)  pd (z) +d (z) p^ (z)  are  strictly 


in  the  unit  cercle 


A-3.2  Regulation  constraints 

A-3.2.1  Passive  regulation 

We  want  to  impose  relation  (18) .  With  the 
hypothesis  on  the  process  and  if  the  sensor  has  a  non 
zero  static  gain,  it  is  necessary  and  sufficient  that  : 

c ( 1 ) *0  (21) 
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So,  we  must  impose  £  factorization  of  the 

compensator  in  : 

c  (z)  -  (z-1 )  m  (z)  (22) 

Moreover,  from  now  on,  we  will  write  the 
sensor  as  k  d(z)  with  : 

d(l)-l,  k*0  (23) 


A-3 .2.2 


impose  : 


Active  regulation 

From  the  preceding  results,  we  must 


b ( 1 ) =0 


(24) 


So,  we  have  the  following  factorization 

b  (z)  =  (z-1 )  n  (z)  (25)  . 


A-3.3  Tracking  constraint 

We  verify  expression  (20)  if  we  impose  identical 
static  gains  for  both  sensor  and  reference.  So  as  in  (23), 
from  now  on,  we  will  write  the  reference  as  k  r(z)  with  : 

r  ( 1 )  *  1  ,  k /0  (26) 


A-3 . 4  Remark 

With  expressions  (22) ,  (23) ,  (25)  and  (26) ,  (6)  roust 

be  rewritten  as 

zm  (z)  e  (z)  =m  (z)  e  (z)  +  (z-1 )  n  (z)  v  (z)  +k  (r  (z)  u  (z)  -d  (z)  s  (z) )  (2  7) 
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So  the  control  is  computed  in  a  recursive  way. 


A-3.5  Stability  constraint 

We  study  the  polynomial  : 

(z-l)o(z)pd  (z)+icd(z)pn(z) 

We  know  already  : 

p_(z)  has  all  its  roots  strictly  in  the  unit  circle 
a 

k,d(l),p  (1)  are  different  from  zero 
n 

-  degree  of  m(z)  is  greater  than  degree  of  d(z) 

-  degree  of  p^fz)  is  greater  than  degree  of  pn(z) 

with  no  more  hypothesis  on  the  process,  we  can  give  a  sufficient 
condition  of  internal  stability  (Proof  in  Appendix  1) . 

If  m(z)  has  all  its  roots  strictly  in  the  unit 
circle,  it  exists  a  vicinity  of  zero  V(0)  such  that  if 
k  is  in  V( 0)-(0),  internal  asymptotic  stability  is  ensured 
if  and  only  if  : 

km ( 1 ) P  ( 1 ) >  0  (stability  condition)  (28) 

with  P(l)  equal  to  the  static  gain  of  the  process. 

Note  that  from  continuity,  the  existence  of  a  vicinity 
of  k  can  be  transposed  on  the  existence  of  a  vicinity  of 
P (z)  as  we  will  see  in  a  next  section,  and  so  this  permits 
the  study  of  robustness  as  it  was  formulated  in  the  introduc¬ 
tion  . 
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A-3.6  Introduction  of  a  non  linearity  on  the  control 


We  will  extend  here  the  results  of  Rouhani  [4] .  we 
introduce  a  non  linear  compensator  defined  as  follows 
(figure  4)  . 


Let  be  the  input  signal  of  the  compensator,  we 

compute  the  control  e^  through  the  expression  : 


N-l 


with  m(z)»  >  m^  z 


(29) 


i =0 
N-i 


fR(x)  a  real  time  varying  function 


V  « 

fn<*> 

*  n  1 

(z-l)m(z) 

7  n  +  1 

FIGURE  4 

-  Non  linear  compensator 

To  study  the  behaviour  of  the  closed-loop  system,  we 
give  an  asymptotic  value  u  to  the  set  point,  we  compute 
a  theoretic  asymptotic  value  e  of  the  control  : 

p  ( 1 )  e'ssu  (30) 

We  suppose  the  disturbances  to  be  bounded  and  the 

—  Np 

processes  to  be  a  M . A .  system  (P(z)=z  p  p^tz)). 
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Then  we  can  say  (Proof  in  Appendix  2)  : 

Let  />  be  the  greatest  modulus  of  the  roots  of 
( z- 1 ) m ( z) zNp+kd ( z) Pn ( z) ,  if  for  any  n  we  have 

|fn(x  +  e)  -e|  <  jr  |x|  ,  k  <  1  (31) 

then  the  non  linear  system  is  stable. 

In  fact  here  (with  the  hypothesis  on  the  disturbances) 
the  stability  is  taken  in  the  sense  of  bounded  input  bounded 
output  (bibo) .  But  if  the  external  signals  (set-point, 
disturbances)  become  constant,  it  will  become  an  asymptotic 
stability  and  verify  relation  (4) . 


A- 4  Sensitivity  of  the  convergence  criterion 


Following  our  introduction  we  are  going  to  study  the 
sensitivity  of  our  preceding  results  to  variations  of  P 
and  Q.  In  fact  given  the  parameters  m,r,k,d,m  of  the 
control  law,  we  are  looking  for  the  set  of  P ,Q  operators 
for  which  the  convergence  criterion  is  satisfied.  To  keep 
the  validity  of  our  approach,  we  will  take  P , Q  in  the 
class  of  linear  time  invariant  processes. 

At  once  let  us  remark  that  only  the  stability  constraint 
uses  hypothesis  on  P  and  Q,  so  we  can  conclude  to 
insensitivity  of  the  tracking  and  regulation  constraints. 

And  from  now  on  we  will  look  at  the  stability  problem. 


A-4.1  Sensitivity  to  Q 

From  expressions  (9)  to  (14)  it  is  easy  to 
conclude  that  for  any  asymptotically  stable  Q,  we  will  have 
internal  stability.  So,  in  fact,  there  is  no  sensitivity  to 
Q  • 


A-4.2  Sensitivity  to  P 

In  the  hypothesis  Hi  we  have  imposed  P  additional 
constraints  to  those  on  Q,  particularly  rational  type  and  non 
zero  static  gain.  The  latter  was  essential  in  regulation 
and  stability  constraints.  So  we  must  impose  variations  of 
P  to  maintain  the  sign  of  the  static  gain.  The  former  was 
a  theoretic  facility  but  it  can  be  relaxed. 


A - 4 . 2 . 1  P  of  rational  type 

In  that  case  we  have  to  find  all  the  pairs  of 
polynomials  (p  (z),  p  (z))  such  that  : 
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G  is  linearly  dependant  on  P. 

Otherwise,  given  the  highest  degree  coefficient 
of  mlz)  p^(z),  from  the  continuity  of  the  coefficients  on 
the  roots,  we  can  say  that  the  set  of  admissible  Gs  which 
represent  polynomials  whose  roots  are  in  the  unit  circle 
is  closed,  bounded  and  connected.  Moreover  from  the  presence 
of  (z-1 ) ,  we  can  say  that  M  is  on  the  frontier  of  this  set. 
Thus  we  have  the  situation  given  by  figure  5. 
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So  from  the  knowledge  of  the  set  of  admissible 
Gs ,  we  can  find  the  set  of  admissible  Ps.  The  first  set  has 
been  studied  by  Markov  ^6}  in  the  continuous  case.  Particu¬ 
larly  we  can't  make  sure  of  convexity  of  the  set,  so  from 
the  linearity  we  don't  know  if  the  set  of  admissible  Ps  is 
connected . 

This  approach  gives  the  roles  of  m,  d  or  k  : 
m  corresponds  to  a  translation,  d  is  very  similar  to  a 
rotation  and  k  to  a  linear  displacement.  Moreover  we  can 
see  the  coupling  between  vicinities  of  k  and  P. 
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A- 4. 2. 2  P  in  the  vicinity  of  a  rational  Pg 

Suppose  the  parameters  of  the  control  law  to 
be  fitted  to  a  nominal  process  Pg  of  rational  type.  We 
are  looking  for  variations  AP  around  Pg,  such  that  we 
have  internal  stability. 

If  we  suppose  P(z)  to  be  an  analytic  function 
outside  a  domain  strictly  contained  in  the  unit  circle  and 
if  we  note  : 

g (z) = (z- 1 ) m (z) +kd ( z) PQ (z)  (3' 

h (z) =kd (z)  A  P  (z)  (3! 

Then  g  and  h  (■—)  are  analytic  in  and  on  the 
unit  circle  and  with  the  Roue he  Theorem  [?]  ,  we  can  say 

that  for  any  process  P  ( z )  =P g  ( z )  *-AP  (  z )  such  that  : 


|g  (elS  )|  >  |h  (eiS  )j 


we  will  have  internal  stability. 


£  [-  T1  ,  n) 


We  have  in  fact  here  another  presentation  of 
the  result  of  Doyle  [8]  in  the  SISO  case. 

A-4.2.3  Application  to  a  polynomial  variation 

Let  us  take  A.P(z)  of  the  form  : 

N 

&P(z)  =  ^  dpjZM_:i  (37) 

j  =0 

Expression  (36)  means  : 

|(eU-l)m(eiB)+kd  (eie)Pg  (ei9)|>  k  Id  (e1 9  )  |  |ap  ( e1® )  (38)| 
V©  £  [  n,  -nj 


But  we  have  with  appendix  3  : 

2  N 

I  A  P  { e  )  |  -  ^  Apj  APi  cos(j-l)© 

3,1-0 

I  i 2  wr-  .  _ 2tf n+  i ,  lisinjN+ne 


i  19,i  a —  2..N+1  .  li  Sin  (N+l 

laPle  >|  APjU-r->?l  — „  T 

3=0 


So  we  can  get  an  upper  bound  of  the  modulus  : 

(  y~  ^p2)2  <■  min  ^e1- 1 )  m  (e1®  ) +kd  (ei8  )  P^e1®^ 

1-0  3  L-*  •+J‘l  k 


Practically  an  FFT  algorithm  will  provide  all 


these  spectra. 


A-4.2.4  Convergence  criterion  sensitivity  index  (CCI) 


Given  a  process  P^  and  the  parameters  of  the 
control  law  we  define  an  absolute  index  by  : 


CCI(P0,m,k,d)= 


lin  I  (ei9 -l)m(eiB)  .  „  ,_i®. 


From  expression  (36)  this  index  gives  an  upper 
bound  on  the  possible  spectrum  variation  to  verify  convergence 
criterion.  So  we  call  it  a  convergence  criterion  sensitivity 
index.  To  insure  robustness,  it  has  to  be  compared  with  an 
equivalent  approximation  index  given  by  the  Pq  model 
estimation  phase. 
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A-S  Summar 


In  this  first  part,  we  came  back  on  the  problem  of 
linear  control.  The  most  important  results  have  been 
reformulated  in  a  very  general  way  :  results  on  the  structure 
of  the  control  law,  results  on  stability  in  the  linear  case 
and  in  a  simple  non  linear  case  and  at  last  results  on  a 
measure  of  the  sensitivity  of  stability. 


B.  THE  SINGLE  INPUT-SINGLE  OUTPUT  (SISO) 

ADAPTED  MODEL  ALGORITHM  CONTROL  (AMAC) 

We  have  just  presented  a  linear  time  invariant  control  law  in  a 
general  fashion.  It  is  an  abstract  approach  which  serves  only  to  ensure 
the  convergence  criterion.  In  an  attempt  to  get  behavior  criteria,  we  are 
going  to  give  a  physical  presentation  through  a  SISO  control  based  on 
the  mathematical  representation  of  the  process  of  paragraph  A-2.1.3: 

s(z)  =  P(z)e(z)  +  Q(z)v(z)  +  w(z)  (1) 

and  the  use  of  adapted  models  of  the  operators  P,Q. 

B-l .  General  SISO  AMAC  Presentation 
B- 1 . 1  Definition  of  the  St rategy 

At  time  n,  given  the  past  measurable  signals,  the  SISO  AMAC  computes 
a  control  such  that  a  predicted  output  of  the  process  is  identical  to  a 
predicted  set  point. 

Taking  the  notation  of  Box  and  Jenkins  [9],  we  write  this: 

snO)  =  un(l)  (2) 

the  prediction  being  here  of  one  point  ahead.  From  the  representation 
of  the  process  (1)  we  decompose  the  predicted  output  into  two  parts:  a 
deterministic  part  which  functionally  depends  on  the  inputs  and  a  non- 
deterministic  part  sun(l)  resulting  from  the  disturbances.  Let  M(en(l), 
e^;  lin)  be  a  model  of  the  operator  P  which  defines  the  deterministic 
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output  from  the  future  and  past  controls,  we  deduce  from  (2)  the  control 
law  as: 


M(en(1),  e£;  Jt  <.  n)  -  un(l)  -  sun(l)  (3) 

Then  to  compute  the  control  en(l)  we  have  three  different  problems: 
inversion  of  M,  estimation  and  prediction  of  the  non-deterministic  out¬ 
put,  and  prediction  of  the  set  point. 

B-1.2  Inversion  of  the  Model  M 

From  its  definition,  M  is  a  model  of  the  process.  Note  that  to 
compute  e^Cl) ,  we  use  this  model  in  a  reversed  way  compared  to  the  physical 
transfer,  so  we  require  M  to  be  invertible  in  the  sense  defined  by  Box 
and  Jenkins  [9]  and  we  call  it  a  deconvolution  model.  Thus  with  the 
linear  time  invariant  hypothesis  the  model  of  the  process  is  taken  linear, 
time  invariant,  asymptotically  stable,  of  rational  type  and  invertible. 

Let  md^  or  Md(z)  be  the  impulse  response  and  the  rational  z-transform 
of  this  deconvolution  model.  We  obtain  from  (3) 


md0,en(1)  “  -  *  "Wl-i  +  Un(1)  "  SUnU) 

i»l 


(4-) 


and  md 


0 


must  be  different  from  zero. 


B-l . 3  Estimation  and  Prediction  of  the  Non-Deterministlc  Output 

From  expression  (1)  the  non-deterministic  output  is  the  sum  of  both 

a  filtered  measured  disturbance  v  and  an  unmeasured  disturbance  w  . 

n  n 

Suppose  we  have  an  estimation  w  of  w  and  a  convolution  model  of  the 

n  n 


1 


measured  disturbance  filter  Q  :  nct(Nc(z)),  then  if  vn(l)  and  vn(l)  are 
predictions  of  measured  and  unmeasured  disturbances,  we  compute: 


sun(D  -  nc0vn(l) 


i*l 


nci‘Vn-H-i  +  °n(1) 


(5) 


So  we  first  need  the  estimation  wn<D  of  the  unmeasured  disturbance  and 
secondly  measured  and  unmeasured  disturbance  predictors* 

We  already  introduced  a  convolution  model  Nc  of  Q.  Let  us  take  also 
a  new  model  mc^McCz))  of  the  process  P.  This  time  we  need  a  model  to  be 
used  in  the  same  way  as  the  process  so  Mc(z)  is  a  convolution  model  com¬ 
pared  with  Md(z),  a  deconvolution  model.  Similarly  to  expression  (1),  we 
compute  the  estimation  w^  by 


a 

n 


s  - 
n 


l 

1*0 


me 


i  n-i 


OO 

l 

i=0 


nc .  • 

l 


n-i 


Now  from  the  past  vn  and  w^,  we  want  to  predict  sun(l).  From  discrete 
parameter  prediction  theory  [10],  v^(l)  ant*  can  comPute<*  w*t*1 

prediction  filters.  Using  3-transforms  they  may  be  expressed  as 

wn(l)(z)  *  Fw(z)  w(z)  - 

(7) 

v(l)(z)  -  Fy(z)  v(z)  *  v(z) 

where  fwn(z),  fwd(z),  fvn(z),  fvd(z)  are  polynomials  in  z,  the  degree  of 
fwd(resp  fvd)  being  greater  than  the  degree  of  fwn(resp  fvn).  Moreover, 
to  be  able  to  predict  the  continuous  component  of  the  disturbances,  we 
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1 

1 


A 


Impose  unit  static  gain  predictors. 

Thus  we  get  the  z-transform  of  sun(l): 

su(l)(z)  =  Nc ( z)  v(z)  +  Fw( z)  w(z)  (8) 

with  Nc'(z)  (ncp  computed  from  Nc  and  through  the  relation: 


nC0  Vn(1) 


QD 

+  l 

i=l 


nc . 
1 


n+l-i 


nc , 


i=0 


n-i 


(9) 


Now  from  z-transform  of  (6)  we  have  the  final  relation: 

su(l)(z)  =  Fw(z)(s(z)  -Mc(z)e(z))  +  (N^(z)  -  Fw(z)Nc(z))v(z)  (10) 


or  equivalently  in  the  time  domain: 

sun(l)  =  fw1*(si-mc1*ei)  4-  (nc!^  -  fwi*nc1)*vi  (11) 

where  *  represents  the  discrete  convolution  operator. 

B-1.3  Set  Point  Prediction 

To  get  a  better  behavior  of  the  closed-loop  system,  at  time  n  we 
need  future  set  points.  In  some  cases  they  are  available  (particularly 
when  there  is  a  hierarchical  control).  But  generally  we  need  a  predictor. 
Let  us  take  it  in  the  form  of  a  rational  filter  Fu(z)  with  fui  as  impulse 
response  and  with  unit  static  gain. 

u(l)U)  -  Fu(z)u( z)  -  u(z)  (12) 

or 


222 


U  (1)  *  fu  *u. 
n  i  i 


(13) 


with  fun(z),  fud(z)  polynomials  in  z,  with  the  degree  of  fud(z)  greater  than 
or  equal  to  the  degree  of  fun(z) . 

B- 1 . 4  Expression  and  Properties  of  the  SI  SO  AMAC  Law 

From  expressions  (3),  (11)  and  (13)  we  get  the  SISO  AMAC  law 

e  (1)  =  — I (fw  *mc .  -  md . )*e  + fu  *u .  -  fw  *s  -  (nc !  -  fw  *nc ,)*v  ] 
n  md^  li  li  li  ii  l  ill 

(14) 


This  prediction  is  used  as  the  future  control  Thus  we  get  the 

z-transform  representation  of  the  SISO  AMAC  law: 

(z*Md(z)  -  FV(z) •Mc(z))o(z)  *  Fu(z) *u(z)  -  Fv(z)s(z) 

-  (Nc'(z)  -  Fw(z)Nc(z) )v(z)  (15) 


Then  we  find  the  expression  of  the  four  polynomials  of  our  general  linear 
time  invariant  control  law: 


c(z)  =  z*Md(z)  -  Fw(z)Mc(z) 
r (z)  =  Fu(z) 
d(z)  =  Fw(z) 

b(z)  *  Nc'(z)  -Fw(z)Nc(z) 


(16) 


Thus  we  can  give  a  physical  interpretation  to  these  polynomials.  Moreover, 
we  see  that  from  a  physical  point  of  view  c,  r,  d  and  b  are  not  mutually 
independent,  but  Md  or  Me,  Nc,  Fu,  Fv  and  Fw  are. 
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The  regulation  and  tracking  constraints  (see  Part  A)  are  verified 
since  we  have  imposed: 

Md(l)  «  Mc(l) 

Fu(l)  =  Fv(l)  -  Fw(l)  -  1  (17) 

Nc'(l)  =  Nc(l) 

The  stability  constraint  (Appendix  4)  can  be  verified  by  a  modification 
of  the  dynamic  of  the  unmeasured  disturbance  predictor  if:  the  different 
models  and  predictors  are  stable,  the  numerator  of  the  deconvolution 
model  has  all  its  roots  strictly  in  the  unit  circle,  the  following  in¬ 
equality  is  satisfied:  Md(l) *P(1)  > 0.  (18) 

Now  assuming  a  perfect  knowledge  of  the  process  and  the  measured 
disturbance  filter: 


Mc(z)  »  P(z);  Nc(z)  *  Q(z) 


(19) 


we  can  write  the  expected  closed-loop  tracking  and  regulation  transfer: 


Sfl(z) 


Fu(z)P(z) 
z  Md(z) 


(20) 


Srv<‘> 

Srv<2> 


Fw(z)Mc(z) 
z  Md(z) 


Nc1 (z)Mc(z) 
z  Md (z 


(21) 

(22) 


So  the  closed-loop  tracking  transfer  is  the  product  of  the  set  point  pre¬ 
dictor  and  the  deconvolution  model  mismatch  of  the  process.  Similarly  we 
get  the  closed-loop  regulation  transfer.  Thus  in  a  perfect  matching, 
the  various  predictors  specify  the  tracking  and  regulation  closed-loop 


transfers. 


The  block  representation  of  the  SISO  AMAC  is  given  by  Figure  6. 


Figure  6.  SISO  AMAC  Representation 

Me  convolution  model  of  the  process 
Md  deconvolution  model  of  the  process 
P  deterministic  part  of  the  plant 

Q  stochastic  part  of  the  plant  (disturbance  process) 

Nc  convolution  model  of  the  disturbance  process 

Nc'  predictive  convolution  model  of  the  disburbance  process 

Fu  set  point  predictor 

Fw  unmeasured  disturbances  (w  )  predictor 


B-2.  SISO  AMAC  Examples 

Following  our  presentation  we  will  present  two  classical  control 


systems  used  whenever  the  convolution  and  deconvolution  model  can  be  iden¬ 


tical. 


with  Md(z)  supposed  to  have  all  its  roots  strictly  in  the  unit  circle. 

Then  if  we  take  a  unit  gain  element  as  a  set  point  predictor,  and  a 
k-step-predictor  zkH(  z)  for  the  disturbance,  we  obtain  the  optimal  control 
system  of  Phillipson  (Figure  7)  which  is  an  improvement  over  the  Smith 
controller. 


w 


Figure  7.  Optimum  Control  System  of  Phillipson 


As  mentioned  by  Phillipson,  this  system  used  in  regulation  is  equivalent 
to  the  Box-Jenkins-Astrom  minimum-variance  control  [9]  or  to  the  Kalman 
linear  regulator  [12]. 

Thus,  this  system  is  essentially  made  for  regulation.  Moreover,  the 
use  of  the  inverse  model  as  controller  since  this  will  be  a  high-pass  filter, 
amplifies  noise,  causes  violent  changes  in  the  control  signal  and  perhaps 
frequent  saturation.  That  is  why  the  AMAC  uses  here  an  adapted  model  and 
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avoids  rapid  changes  in  set  point  thanks  to  the  set  point  predictor. 
Predictor  H  can  be  easily  computed  when  the  disturbance  can  be  described 
as  the  output  of  a  known  rational  filter  whose  input  is  an  independent 
zero-mean  random  sequence.  But  to  verify  internal  stability  we  must 
not  forget  the  constraints  on  the  denominator  of  H.  Here  Phillipson 
suggests  the  use  of  exponential  smoothing  for  prediction  to  solve  the 
problem.  That  way,  we  can  answer  satisfactorily  the  output  regulation 
but  not  so  properly  the  output  tracking.  The  model  predictive  heuristic 
control  which  follows  attempts  to  answer  the  two  questions  introducing 
a  set  point  predictor  and  deducing  the  disturbance  predictor. 

B-2 . 2  Model  Predictive  Heuristic  Control  (MPHC)  [13] 

We  give  here  a  simplified  study  of  this  method;  the  very  general 
study  can  be  found  in  [14].  Suppose  no  measured  disturbance  (MPHC  can 
be  extended  to  this  case)  and  a  convolution  model  having  a  moving  average 
(MA)  structure  with  all  its  roots  strictly  in  the  unit  circle,  we  take: 


M(z)  =  Me  ( z)  =  Md(z) 


For  the  set  point  and  disturbance  predictors,  we  choose 


F  (z) 
u ' 


l-G(l) 

1  -  z  ^G(z) 


F 

w 


(z) 


l-G(z) 

1  -  z  ^"G(z) 


(24) 


(25) 


where  C(z)  is  a  nonzero  static  gain  transfer  such  that  Fu(z)  and  Fw(z) 
satisfy  stability  conditions. 
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Then  from  (15)  we  get  the  MTHC  law: 


)  *M(z)  *e(z)  =  - -  [  (1  -  G(l)  )u(z)  -  (1  -  G(z))  s(z)  ] 

1-z  G(z)  1-z  G(z) 

(25) 


(z-l)*M(z)*e(z)  +  s ( z)  =  (1-G(l))u(z)  +  G(z)s(z) 


(26) 


Let  us  develop  the  strategy  of  this  relation. 

Both  terms  are  similar  to  outputs.  We  call  the  left-hand  term  a 
predicted  output  Sp(z)  and  the  right-hand  term  a  reference  output  s^(z). 
From 


Sp ( z)  =  zM(z)e(z)  +  (s(z)  — M(z)e(z))  (27) 

we  define  the  predicted  output  as  the  output  of  the  model  at  time  (n+1) 
corrected  of  the  estimation  w(n)  of  the  disturbance 

sp(n+l)  =  sM(n+l)  +  w(n)  (28) 

with  s  (n+1)  output  of  the  model  with  a  predicted  input  e  (1).  The  ref- 
M  n 

erence  output  sR(z)  is  given  by  a  trajectory  connecting  the  past  outputs 
of  the  process  to  the  present  set  point. 

sR(n+l)  =  (1  -  l  gi)un  +  l  sn_.  (29) 

i“0  i=0 

Thus  the  MPHC  strategy  consists  in  computing  future  inputs  such  that 
predicted  outputs  are  on  a  connecting  trajectory.  Its  block  representation 
Is  given  by  Figure  8. 
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l-C(l) 

+ 

1 

1-z  ^G(z) 

M(z) 

Figure  8.  MPHC  Representation  with  G(z)  as  Connecting  Trajectory  Generator 


From  part  A,  we  will  satisfy  the  convergence  criterion  if  M(z)  has  all 
its  roots  strictly  in  the  unit  circle  and  (l-G(l))  is  taken  as  the  stability 
coefficient.  But  again,  the  transfers  are  not  independent,  for  instance  in  a 
perfect  modeling  we  have: 


sa(z) 


l-G(l) 
z  -  G(z) 


S  (z)  =  1  - 
rw  z  -  G(z) 


z-1 

z-G(z) 


l-G(l)  Sa(z) 


The  regulation  transfer  is  the  discrete  differentiation  of  the  tracking 
transfer.  Moreover,  if  the  model  does  not  verify  the  stability  condition, 
the  strategy  must  be  seriously  questioned  but  it  has  been  extended  to  this 
case  by  the  introduction  of  the  notion  of  adapted  model  [14]. 


B-3.  Choice  of  the  Deconvolution  Model 


We  have  seen  that  the  most  general  linear  time  invariant  control  law 
contains  five  independent  physical  components.  Theoretically  each  can  be 
obtained  from  a  modeling  (system  or  spectrum).  But  the  deconvolution 
model  is  a  special  case  because  its  use  is  not  a  physical  one.  We  are 
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going  to  show  where  the  problem  is  and  how  to  solve  it. 

B-3.1  Terms  of  the  Problem  and  Mathematical  Solution 
Let  Md(z)  be  this  deconvolution  model 


_  mdn(z) 

Md(z)  “  ^SddUT  (32) 

with  mdn(z),  mdd(z)  polynomials  in  z  such  that  with  expectation  (4)  degree 
of  mdd(z)  is  equal  to  degree  of  mdn(z). 

Mc(z)  is  the  knowledge  of  the  process,  i.e.,  the  convolution  model: 


Mc(z) 


men( z) 
mcd(z) 


(33) 


with  mcn(z),  racd(z)  polynomials  in  z,  with  degree  of  mcd(z)  greater  than 
degree  of  mcn(z).  The  differences  between  these  models  are  in  their  use. 
Let  e,s  be  the  input  and  the  output,  we  write 


s(z) 


mcn(z) 

racd(z) 


e(z) 


(34) 


similarly  to  the  process,  but: 


e(z) 


mdd(z) 
md  n  (  z ) 


s(z) 


(35) 


is  a  reversed  relation  compared  with  the  process. 

As  we  want  a  stationary  control  law,  following  Box  and  Jenkins  [9], 

we  have  to  improve  the  stability  of  both  transfers  — C-^— c  and  . 

mcd(z)  mcn(z) 

The  former  can  be  ensured  from  the  stability  of  the  process.  But  the 
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latter  has  no  physical  significance  and  we  have  seen  that  we  must  impose 
mcn(z)  to  have  all  Its  roots  strictly  in  the  unit  circle. 

As  mdd(z)  has  no  constraint  in  the  deconvolution  use,  we  can  take: 

md(z)  =  mdd(z)  =  mcd(z)  (36) 

Moreover,  to  get  a  zero  static  gain  compensator,  with  expression  (16), 
we  must  impose: 

mdn(l)  =  mcn(l)  (37) 

Thus  the  problem  is:  knowing  the  model  of  the  process  mcn(z),  how  to 
choose  mdn(z)  such  that  it  keeps  the  significance  of  a  model  and  it  satisfies 
the  stability  condition. 

If  mcn(z)  has  all  its  roots  strictly  in  the  unit  circle,  we  take 
obviously : 


mdn(z)  =  men(z)  (38) 

So  the  real  problem  occurs  when  mcn(z)  has  roots  on  both  sides  of  the 
unit  circle.  Let  us  factorize  mcn(z)  into: 

mcn(z)  3  min( z) *mon(z)  (39) 

where  rain(z)  has  all  its  roots  strictly  inside  the  unit  circle,  mon(z) 
has  all  its  roots  strictly  outside  the  unit  circle.  We  don't  deal  with 
unit  modulus  roots.  As  ndn(z)  is  used  as  a  denominator,  let  us  consider: 


mid(z) 


md  ( z ) 
min( z) 


(40) 
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mod(z) 


1 

mon(z) 


(40) 


mid(z)  is  a  holoraorphic  function  defined  outside  a  domain  strictly  con¬ 
tained  in  the  unit  circle,  so  its  Laurent  expansion  in  the  vicinity  of 
the  unit  circle  is: 


mid(z)  =  £  mid  z  ^  (41) 

j=0  J 

mid(z)  corresponds  to  a  causal  impulse  response  and  so  has  a  physical 
significance. 

Inversely,  mod(z)  is  a  holomorphic  function  defined  in  a  domain 
strictly  containing  the  unit  circle,  so  its  Laurent  expansion  in  the  vi¬ 
cinity  of  the  unit  circle  is 


mod(z)  =  £  mod.zJ 


j=0 


j 


(42) 


Thus,  mod(z)  can  be  considered  as  corresponding  to  a  noncausal  impulse 

response.  And  so  expression  (35)  or  (3)  implies  the  knowledge  of  the 

future  outputs:  e  (1)  is  functionally  dependent  on  u  (k),  su  (k),  keN. 

n  n  n 

Precisely,  from  (3)  we  get: 


IP'md(z)  =  mod(z)  [u(l)  (z)  —  su(l)  (z)  ) 


(43) 


e  (1)  depends  on  the  term 
n 


I  rood  (u  (j)  -su  (j)) 
j=0  3 


(44) 
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i.e.,  for  the  one  step  ahead  prediction  input,  we  need  j-step  predictions 
of  the  set  point  and  the  non-deterministic  output  for  all  positive 
integers  j.  This  is  a  mathematical  result;  the  physical  problem  is  that 
predictions  are  not  real  values.  Thus  the  strategy  of  the  SISO  AMAC  cannot 
be  totally  ensured. 


We  are  going  to  show  how  this  mathematical  solution  can  be  used  to 


design  the  control  law. 

To  simplify  the  statement,  we  will  suppose  no  measured  disturbance 
and  as  suggested  by  Phillipson  and  MPHC  applications,  we  take  exponential 
smoothing  for  prediction: 


Fu(z) 


F  (z) 
w 


1  -  tz 


1  -  rz 


with  t,r  called  tracking  or  regulation  coefficients.  Moreover,  as  there 
is  no  problem  on  the  model's  denominator,  we  suppose  an  MA  model  (with 
p  the  number  of  coefficients) 


Mc(z)  =  ECU! 

_P 


with  mc(z)  factorized  in  mi(z)*rao(z). 

So  we  will  work  with  the  block  representation  given  by  Figure  9  and 
the  AMAC  law  given  hy  the  relations 


"V„(1)  ■ 


-  I  "Vn+l-l  +  ”„<»  '  ‘“.W 


(47) 


un(j)  “  un(l)  -  tun_1(D  +  d-r)un;  u0(l)  -  uQ 
sun(j)  =  sun(l)  =*  t  su^d)  +  (l-r)wn;  sQ(l)  =  wQ 


Figure  9.  AMAC  Representation  for  Study 


B-3.2  Direct  Application  of  the  Mathematical  Solution 

From  the  factorization  of  mc(z),  let  mi^  be  the  impulse  response 
corresponding  to  the  roots  inside  the  circle  and  mod^  be  the  noncausal 
impulse  response  corresponding  to  the  inverse  of  the  factor  containing 
the  roots  outside  the  unit  circle.  We  write  the  control  law  as: 


min*)en(l)  “  "  [  minjen+1-j  +  ^  modj(un(j)  ‘  sun(j)) 


But  with  our  constant  prediction  we  have: 

min*en(l)  -  -  l  min*en+1_j  +  (  £  mod^ ) (ur(1)  -  sun(l) )  (49) 

With  exponential  smoothing  for  predictions  the  mathematical  solution  gives 
a  deconvolution  model  which  has  among  the  roots  of  the  convolution  model 
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only  those  which  are  strictly  in  the  unit  circle. 


From  (20),  (21) 


sa(z) 


1-t  mon( z) 

z-t  0  ... 

z  mon(l) 


S 

rw 


1 


1-r  mon(z) 

z-r  0 

z  mon(l) 


(51) 


With  0  the  number  of  roots  of  mon(z).  Thus,  the  prescribed  behaviors  can 
be  followed  only  after  the  response  time  of  mon(z). 


B-3.3  k-Step  Prediction 

Our  mathematical  presentation  tells  us  that  to  compute  en(l)  we  need 
further  predictions.  So  one  idea  is  to  rewrite  the  AMAC  strategy  as 

s  (k)  -  u  (k)  (52) 

n  n 


with  no  a  priori  distinction  between  the  models. 
Then  (47)  gives 


1  me.  e  (k-i)  -  -  [me.  e  +  un(l)  -  su^l)  (53) 

i=0  i=k 


Thus  the  control  e  (1)  depends  on  the  predicted  inputs  e  (j)  and  we  have 
n  n 

to  solve  a  linear  system  with  k  unknown  quantities.  To  get  a  unique  solu¬ 
tion  one  could  introduce  a  criterion  on  the  predicted  inputs. 

Let  us  look  for  a  solution  linearly  dependent  on  the  right  term  as: 


fe  (k)l 
n 

• 

en(i) 

n 

k 

P 

T  me . • e  .  . 

. L,  i  n+k-i 
i“k 


+  u  (l)-su  (1)) 

n  n 


(54) 
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but  the  deconvolution  model  can  be  easily  computed. 

This  k-step  prediction  strategy  can  be  extended.  Given  a  set  I  of 
positive  inters,  we  impose: 

s  (k)  =  un(k)  ,  V  ke  I  (60) 

This  leads  as  previously  to  a  linear  system  whose  unknown  quantities  are 
the  predicted  inputs.  It  can  be  solved  in  various  ways,  but  the  solution 
must  give  a  deconvolution  model  satisfying  the  stability  conditions  [14]. 

The  advantage  of  such  an  approach  is  in  the  constrained  control 
case:  let  ft  be  the  time  invariant  set  of  admissible  inputs,  we  write  the 
extended  k-step  prediction  strategy  as: 

Min  J(s  ( k)  —  u  (k) ;  kel)  (61) 

en(j)efi  n 

where  J  is  a  criterion. 

This  optimization  problem  gives  predicted  inputs  satisfying  the  con¬ 
straints.  Thus,  one  can  expect  to  get  a  better  behavior  owing  to  the  fact 
that  predicted  constraints  are  taken  into  account. 

B-3.4  Choice  from  Behavior  Analysis 

The  deconvolution  model  defines  the  tracking  behavior  with  the  following 
transfer  obtained  in  a  perfect  modeling  hypothesis. 


S.<*> 


Fu(z)P(z) 


z  Md(z) 


(62) 


When  we  can  take  Md(z)  equal  to  P(z)  the  set  point  predictor  plays  the 


same  role  as  the  reference  model  in  the  MPHC,  so  we  can  extend  here  the 


ideas  of  Rouhani  [4]. 

B-3.4  Pole  Placement 

One  can  Impose  direct  pole  placement.  In  that  case,  from  the  speci¬ 
fied  polynomial  pp(z)  we  get  the  deconvolution  model  as: 


M d(z) 


min(z) pp(z) 
rad(z) 


(63) 


because,  in  perfect  modeling,  we  have 


S  (z)  =  Fu<z)  mon<z) 
a'  z  pp(z) 


(64) 


This  method  is  very  simple  when  the  factorization  (39)  is  known.  If  not, 
this  problem  may  be  very  difficult  to  solve  numerically  in  particular 
when  there  Is  a  great  number  of  roots. 


B-3.5  Optimization  Criterion 

Another  natural  criterion  is  the  minimization  of  a  quadratic  distance 
between  the  expected  and  the  actual  responses  to  a  set  point  sequence: 


J(Md(z)) 


+TT 


Fu(e16) 

10 


M  (e^) 

-s-le-)  u(ei0> 
Md(e16) 


2 

d0 


(65) 


^6 

where  u(e  )  is  a  specified  function. 

This  is  equivalent  to  a  distance  between  deconvolution  and  convolu¬ 
tion  models.  Thus,  if  we  write 
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Md(z) 


(66) 


mdn(z) 

md(z) 


Mc(z) 


roc.n(z) 

rad(z) 


The  problem  is  to  approximate  the  polynomial  mcn(z)  whose  roots  are  on 
both  sides  of  the  unit  circle  by  a  polynomial  mdn(z)  whose  roots  are 
inside  the  unit  circle  and  (65)  can  be  rewritten  as 


J(mdn(z) ) 


•Hr 


1  II 

» 

,  ie, 

.  mcn(e  ) 

-7T 

mdn(ei0) 

dF(0) 


(67) 


with  dF(0)  a  positive  measure  and: 

mdn(l)  =  racn(l)  (68) 

Such  a  criterion  and  constraints  can  be  computed  by  the  Jury-Astrom  al¬ 
gorithm  [15]. 

This  method  gives  a  deconvolution  model  which  depends  only  on  the 
convolution  model  Me  and  the  set  point  predictor  (tracking  coefficient  in 
the  exponential  smoothing  case).  Those  computations  may  be  numerically 
easier  than  polynomial  factorization,  but  have  to  be  done  again  if  the 
predictor  is  changed. 

B-3.6  Cone 1 us ion 

The  deconvolution  problem  can  be  solved  using  a  prediction  strategy. 
This  leads  to  a  simple  method  but  not  manageable  results  in  the  k  step 
prediction  case  or  to  more  difficult  computations  as  polynomial  factori¬ 
zation  or  constrained  nonlinear  optimization  if  we  want  to  have  more 
manageable  results. 
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B-4.  Summer 


We  have  shown  that  a  very  general  implementation  of  a  linear  time 
invariant  control  law  can  be  done  by  the  adapted  model  algorithm  control. 
This  method  uses  five  independent  physical  entities:  two  non-deterministic 
signal  predictors  which  can  be  deduced  from  disturbance  modelization; 
two  mathematical  representations  of  the  process  behavior  which  are  also 
given  by  modelization;  a  set  point  predictor  which  can  be  deduced  from  the 
control  law  specification.  The  problem  is  complicated  by  the  fact  that 
one  of  the  representations  is  used  in  a  nonphysical  way  and  thus  has  to 
be  adapted  by  a  further  prediction  strategy. 
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APPENDIX  1 


Proof  of  the  linear  stability  theorem . 

Let  C  (z) ,A(z) ,B (z)  be  three  polynomials  with  the 
relation  : 

C.  ( z)  =  (  z- 1 )  A  (  z) +k  B(z)  (A-l.l) 

with  degree  of  A(z)  greater  than  degree  of  B(z) 

„ B ( 1 ) /0  (A-1.2) 

.  the  highest  degree  coefficient  of  A(z) 

_  A ( 2 )  has  all  its  roots  strictly  in  the  unit  circle 
-  A ( z ) ,  B(z)  with  real  coefficients. 

We  will  show  that  it  exists  a  vicinity  of  zero  V(0)  such 
that  if  k  is  in  V 10) -  (0),  the  roots  of  Ck(z)  are  strictly 
in  the  unit  circle  if  and  only  if 

k  aQ  B(l)>  0  (A-l .3) 

Proof  :  we  use  continuity  results  :  the  roots  of  a 
polynomial  of  the  complex  variable  and  the  maximum  of  their 
moduli  are  continuous  functions  of  its  coefficients,  the 
highest  degree  coefficient  being  taken  different  from  O. 

So  we  have  : 

1  -  For  any  real  k,  Ck<z)  has  (da+l)  roots  if  is 
the  degree  of  A(z). 

This  is  Appendix  1  of  the  report  which  comprises  Appendix  A  of  the  whole  report 
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2 


-  For  k  equal  to  zero,  the  roots  of  Cq(z)  are  : 
the  d  roots  of  A(z)  which  are  strictly  in  the 
unit  circle 

The  simple  root  :  equal  to  1 . 

3  -  From  the  preceding  continuity  properties,  it 

exists  a  vicinity  V{0)  of  zero  such  that  for 
any  k  in  this  vicinity,  C(z)  has  d  roots 
strictly  in  the  unit  circle. 

4  -  Let  us  study  the  last  root. 


If  the  modulus  is  greater  or  equal  to  one,  the  root  is 
real  because  it  is  alone  outside  the  unit  circle  and  the 
coefficients  of  the  polynomial  C(z)  are  real.  So  let  us 
consider  the  polynomial  C(z)  of  the  real  variable,  its 
only  root  greater  than  one  exists  iff  : 

C(  1  )  C(x)  >  0  ( A- 1 .4) 

for  large  x  greater  than  one. 

But  here  we  have  : 

C(  1 )  =k  B  ( 1 )  (A- 1  .5) 

and  the  sign  of  C(x)  for  large  x  is  the  one  of  a q  so 
there  is  no  root  if  : 

kB ( 1 )  aQ  >  0  (A-l .6) 

Remark  :  From  the  hypothesis  on  A(z)  ,  the  signs  of 
A ( 1 )  and  a^  are  the  same. 
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APPENDIX  2 


Proof  of  the  non  linear  stability  theorem. 

Let  us  take  the  notations  : 

Etn)*  -  (  (e  -e)  .  .  .  (e  -e)  .  .  .  (e  -e)  )t 
n-M  n-N  n 

m*"  =  (0  .  .  .  0  .  .  .  m. )  t 

N  U 

Stn)1  =  ( ( s  -u) . . .  (s  -u) )  t 
n-N .  n 
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with 


*n-N, 


<3. 


ln-N, 


(1 


1) 


(0 


01) 


t.  rp  t 


With  :  H  =  identity  matrix 


M  =  N  +N 
P  d 


L  =  N  +N 
q  d 


N  number  of  coefficients  of  m(z) 


N  number  of  coefficients  of  d(z) 
d 


N  number  of  coefficients  of  r(z) 
r 


N.  number  of  coefficients  of  b(z) 
b 


those  notations,  we  write-the  compensator  relation  (29) 


e  *  =  f  (-7-  f  y  +mt(n  -  IT)  (E(n)+e  t  )]  ) 
n+ 1  n  m„  1 ' n  0 


(A-2.2) 


-the  compensator  input  (27)  : 

yn*k[rt  (U  (n)  +u  tQ) -dfc  (S  (n)  +u  tQ)j-btV(n)  (A-2.3) 
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-and  the  process  output  (2)  : 


S(n)«  (P  E(n)  +  C  V(n)+W(n)  (A-2.4). 

From  the  equality  of  static  gains  of  the  sensor  and 
the  reference,  we  have  : 

u  rtt^=u  d^Q  (A-2.5) 

So  : 

y  =  k  rtU (n) -kdfc  (P  E (n) - (kdt  (Q  +bt) V (n) -kdfcW (n)  (A-2.6) 

n 

And  : 

e  =f  (e+”  f(mt(n  -TT) -kdfc  |-p  )E(n)+  krtU(n) 
n+1  n  mo 

~(kdtC?  +bt)  V  (n) -kdfcW  (n)j  )  (A-2,7) 

Thus,  we  have  a  state  representation  of  the  control 

E(n+1)=F  (A  E(n)+x  t,+e  t_)-eT  t_  (A-2.8) 

n  n  1  U  u 

with  : 

(A-2 .9) 


(A-2. 10) 
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irw 


X  «—  (  krtu  (n)  -  (kd^  +bt)  V  ( n )  -  kd tW  (n )  ) 
n  m0 


(h-2 . 11) 


Now  we  remark  that  is  a  companion  matrix  associated 
to  the  polynomial  : 


P  t 

(z-l)m(z)z  +k  d  p  (z) 

n 


So  we  have  a  new  result  of  stability 


Let  f>  (A)  be  the  spectral  radius  of  the  matrix  /^  , 

if  f  satisfies  the  following  inequality  : 
n 


|fn(x+.U.|  < 


x  I  ,  k  <  1 


(A-2 .  12) 


Then  the  system  with  the  non  linearity  f  is  stable. 


Proof  : 

If  the  linear  system  is  asymptotically  stable,  the 
spectral  radius  of  /&  is  less  than  one,  then  it  exists  a 
consistent  norm  of  which  is  less  than  one  ^5]  . 

For  that  norm,  we  have  : 


x  ||  <  (M)  .  ||  x  ||  ,  ( R)  C  1 


II  AiE  (n)  +xRt  1  ||  <  f>(Pi)  l|E  ( n )  W  +  IxJ  .lltjl 


But  from  (42) ,  we  have  also  : 


(A-2. 13) 


(A-2 . 14) 


Fr  (  E  (n)  tx^tj+etQl-etQ^  <~f(  Ai  )  E  (n)  +x^t  (A-2. 15) 
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So  : 


lxnl  '  , 

||  E  (n+ 1 )  (|  4  k(  llE(n)H  + 


(A-2. 16) 


Then  if  is  bounded  i 
and  unmeasured  disturbances  are 
our  proof . 


e.  set-point,  measured 
bounded,  we  can  conclude 


;  4 


P 

> 


Id 

* 

i: 
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APPENDIX  3 


Spectral  analysis  of  the  matrix  (cos(i-j)©  ). 


From  the  relation  : 


cos  (i- j  )e>  =cosie. cosj©  +sini©  .  sinjff  , 

we  can  write  : 

A 


(A-3.1) 


.cos(i-j)p  .|=(cosi©  |(.cosj©  .)+|sini© 


(  . sin j©  . ) 

(A-3 .2) 


So  we  see  that  the  matrix  is  semi  definite  positive 
with  only  two  positive  eigen  values.  We  are  looking  for  the 
eigen  vectors  as  a  linear  combination  : 


sinj  & 


x (cos je  }  + 


We  have  to  find  x  and  y  such  that  : 


• 

M 

2 

M 

• 

cosi  0 

(x  .  y  cos  j  8 

+  ^  cosje  sinj©  ) 

xy 

cosi  e 

u- 

ii 

o 

j=0 

• 

+ 

= 

+ 

• 

M 

M  2 

• 

sini  © 

(x.  V  cos  j© 

sinj©  +  \  sin  j©  ) 

y 

sini© 

• 

j  =0 

j=0 

ft 

(A-3 


We  deduce  the  expressions  : 
xA+B*xy 

xB+C*  y  (A-3. 4) 


3) 
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r 


with  A- ^  COS  jP 

j-o 


B»y  sinj&  cosjO 

j-o 


C-^~  sin 

j=0 


We  get  : 


x  B+x (C-A) -B=0 


A-C+ ( (C-A) 2+4B2)2 


A+C+  (  (C-A) 2  +  4B2) 2 


but  : 


A+C-M+l 


A-C-y  COS2  j®.  * 

/  I  sin  Q 

j-0 

n  -  Is  in  (M+l  )& 


2B 

j-0 


cos  j$  * 1 

A-  1  i  ■  ■  ■ 

.  M  ^ 
sin'j  9 


M- 
COSJ  © 


sinj© 


_M+ 1  1  sin (M+l ) & 
2  -2  sin© 


(A-3.5) 


(A-3 .6) 


(A-3.7) 


(A-3. 8) 


(A-3. 9) 
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APPENDIX  B 


APPLICATION  OF  ROBUSTNESS  ANALYSIS  TO  SIMPLE  MISSILE  EXAMPLE 

The  following  example  applies  the  robustness  analysis  of  Section  III 
to  a  single  input-single  output,  second  order  model  of  missile  longitudinal 
dynamics  (see  Section  VIII).  After  some  preliminary  analysis  to  obtain 
the  plant  and  model  z-transform  used  in  the  robustness  analysis,  we 
analyze  the  stability  of  the  controlled  plant  when  some  parameters  (plant 
parameters,  sample  rate,  reference  rate,  etc.)  are  varied. 

1.  Preliminary  Analysis 

As  shown  in  Section  III,  the  stability  of  the  (sampled)  controlled 
plant  is  determined  by  the  zeros  of  the  expression 

(z-l)H(z)  +  (l-a)H(z)  (B.l) 

where  H(z)  is  the  z-transform  of  the  plant  model,  H(z)  is  the  z-transform 
of  the  plant,  and  a  is  the  rate  of  the  first-order  output  reference  trajec¬ 
tory  (0<a<l).  The  MAC  controlled  system  is  stable  if  the  zeros  of  (B.l) 
lie  within  the  unit  disc.  In  this  section  we  derive  expressions  for  H(z) 
and  H(z)  in  order  to  analyze  the  zeros  of  (B.l)  in  the  next  section. 

Suppose  the  real  plant  is  given  by  a  linear,  time  invariant  state 
description 

x  =  Ax  +  Bu  (B.2) 

y  =  Cx  (B.3) 
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If  this  system  is  relaxed  at  t  =  0  (i.e.,  x(0)  =0)  and  if  the  input  u  is 
constant  over  the  intervals  [kA,  (k+l)A)  of  length  A  for  k=  0,1,2,..., 
then 


y(kA)  =  l  Ce^V^e^-DBu^^j  (B.4) 

0 

In  (B.4)  it  is  assumed  that  all  the  eigenvalues  of  A  have  negative  real 
parts  so  that  the  inverse  A~*  is  well  defined.  The  z-transform  corres¬ 
ponding  to  (B.4)  is 

H(z)  =  l  Ce^V^e^DBz'^  (B.5) 

£= 0 

which  is  equivalent  to 


H(z)  =  C(I  (B.6) 

We  assume  that  the  plant  model  used  by  MAC  has  a  form  similar  to  the 
plant  (B.2),  (B.3)  in  the  following  way.  The  model  uses  a  discrete-time 
impulse  response  computed  from  a  linear,  time-invariant  state  description 
with  matrices  A,  B,  C  corresponding  to  A,  B,  C  in  (B.2),  (B.3).  The  time 
interval  is  the  same  as  above,  namely  A.  However,  the  model  has  the 
following  differences.  The  infinite  sum  in  (B.4)  is  truncated  after  a 
finite  number  of  terms.  In  addition,  A"*(e^-I)  is  approximated  by  A*I. 
With  these  assumptions  we  obtain  the  following  expression  for  H(z). 

H(z)  =  l  Ce^BAz-*  (B.7) 

£=0 
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I 


This  is  equivalent  to 

H(z)  =  C(I  -  z"1e^A)_1(I  -  z-(M*1)e^N+1))B  (B.8) 

In  the  examples  we  will  assume  that  A,  B,  C  have  the  following  form 


(B.9) 

(B.10) 

(B.ll) 

The  expressions  for  A,  B,  C  will  be  the  same,  with  Z  ,  M  ,  M  replacing 

W  01  St] 

Zw’  M  ’  Msq  in  (B-9>  -  (B>11)  • 


In  order  to  compute  H(z)  and  H(z)  it  will  be  necessary  to  compute 


as  given  in  Section  VIII. 


A  = 


'Z„  1 

M  0 

a 


B  = 


'0  ' 


lMsqj 


c  =  (1  0) 


M  =-281.11,  which  gives 

6  =  -.7434,  u  =  12.2220 
The  corresponding  H(z)  is 

H(z)  =  -(1.1807)z(z+ .9490)(z2  - (.6342)z  + (.8620))"1 

which  has  zeros  at  z  =  0  and  z  =  -.9490,  both  stable  as  they  should  be. 

The  model  plant  used  by  MAC  in  this  case  will  have  the  same  para¬ 
meters,  i.e.,  ZW  =  ZW,  Ma  =  Ma,  Msq  =  Msq  and  with  N  =  49,  A=  .1.  The  cor¬ 
responding  H(z)  is 

H(z)  =  -(2.0071)z'49(z2-  (.6342)z + ( .8620) )_1(z50 + (.0280)z-  (.0056)) 
Note  that  z  +  .0280z-  .0056  can  be  wnttui  as 

z50+  .0280(z  -  .2000) 

Since  |.0280|  <  h,  Rouche's  theorem  (see  Ahlfors,  1966)  implies  that 
z^9  +  .0280(z-  .2000)  and  z"  have  the  same  number  of  roots  in  the  unit 
circle.  In  other  words,  all  zeros  are  inside  the  unit  disc;  hence,  all 
zeros  of  H(z)  are  stable. 

The  controlled  system's  stability  depends  on  the  zeros  of 

(z-l)H(z)  +  (l-a)H(z)  =  (l-a)(-1.1807)z(z  +  .9490)(z2  -  .6342z+  .8620)'1 
+  (z-l)z'49(z2  -  .6342z+ .8620)'1(z50+ .0280z-  .0056)(-2 .0071) 

The  zeros  of  this  function  are  the  same  as  the  zeros  of  the  polynomial  f(z) 
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f(z)  *  (l-a)(z+  .9490) (-1.1807)  +  (z-l)(z50  +  .0280z- .0056)(-2.Q071) 


which  is  equivalent  to 


f(z)  »  z51(l-a+1.6999)  +  z3U((l-a)( .9490)  -  1.6999)  +  ( .0476) (z-1) (z- .2000) 


.50 


Let  fj(z)  and  f2(z)  be  defined  as 


fj(z)  =  (l-ct+1.6999)z51  +  ((l-aU.9490)  -  1.6999)z50 


f2(z)  *  ( -0476) (z-1) (z- .2000) 


Thus,  f(z)  =  fj(z )  +  f2(z).  For  |z|=l,  it  is  straightforward  to  check  that 
for  a  <  1 


lVz)i <  lfi^)l 

Hence,  by  Rouche's  theorem,  f(z)  has  the  same  number  of  roots  as  fj(z) 
within  the  unit  disc.  It  is  easy  to  see  thatwfj(z)  has  all  its  51  roots 
inside  the  unit  disc  for  0<a<l.  Thus,  the  controlled  system  is  stable 
for  all  0<a<l  in  this  case. 


2.2  Example  2 

In  this  example  the  plant  parameters  are  Z  =-.27877,  M  =-64.928, 

n  O 

M ■  -43.285  which  gives  6  =  -.1394,  ui  =  8.0566.  The  corresponding  H(z)  is 
H(z)  =  -(.0087)z(z-  .9710)(z2- (1.9721)z+  .9725)"1 


The  model  plant  used  by  MAC  in  this  case  will  have  parameters  the 


same  as  for  the  previous  case,  Z  = -1.4868,  M  = -149.93,  =  -281.11  and 

w  a  sq 

with  N  =  49,  &=  .1.  As  before,  R(z)  is  given  by 

fl(z)  =  -(2.0071)z'49(z2-  .6342Z+  .8620)_1(z50  +  .0280z-  .0056) 

The  stability  of  the  controlled  plant  depends  on  the  zeros  of  the  poly¬ 
nomial  f(z)  given  by 

f(z)  =  (-2.0071-  k(.0087))z53 

+  (5.9653  + k(.0140))z52 

+  (-5.9101-  k(.0126))z51 

+  (1.9519-  k(.0073))z50 

+  (- .0562 )z4 

+  ( . 1782 )z3 

+  (- . 1989)z2 

+  ( . 0877 ) z 
+  (-.0108) 

where  k=  1-a  and  a  is  the  parameter  in  the  reference  trajectory,  0<a<l. 
To  test  the  stability  of  the  zeros  of  f(z)  it  is  necessary  to  use  a  num¬ 
erical  test  such  as  the  Schur-Cohn  test  (see  Tretter,  1976,  p.  104)  to 
test  whether  any  zeros  of  f(z)  are  outside  the  unit  disc.  For  example, 
there  are  roots  on  or  outside  the  unit  disc  if  ka0,  1-e"*  (=  .6321),  but 
all  roots  are  inside  the  unit  disc  for  k  =  1. 


APPENDIX  C 


DERIVATION  OF  THE  OPTIMUM  IMPULSE  RESPONSE 
FUNCTION  FOR  NONMINIMUM  PHASE  SYSTEMS 

1 .  Derivation  of  the  Rlccati  Equation  for  p  =  0 
The  optimum  control  problem  of  Section  IV  is: 

C(k+1)  =  UC(k)  +  bju(k)  (C.l) 

J(n)  =  -*^T(n)QC{n) CT( j)QC( j )  (C.2) 

j=0 

N 

Since  dim  [b^,  Ub^,...,  U  bj]  is  equal  to  N+l,  the  system  (C.l)  is  con¬ 
trollable. 

The  cost  function  J(n)  has  the  peculiarity  of  not  depending  explicitly 
on  the  input  u'(j).  But  in  order  to  be  able  to  use  the  common  discrete 
version  of  LO,  one  must  require  that  the  cost  associated  with  the  input 
u'(j)  be  strictly  positive  (Sage  and  White,  1977;  Murata,  1977).  That  is, 
the  cost  function  must  be  of  the  type: 


Jr(n)  =  «T(n)CK(n)+>s  l  {?'  (j)QS(j)  +  pu,d(j)},  with  p>0 
c  j=0 


n-1 


-T 


,  1 2  /  4 1 


In  practice  it  is  always  possible  to  "correct"  J  with  an  additional 


n-1 


o 

term  p  j  u1  (j),  where  p  is  a  small  positive  number,  and  obtain  Jc-  By 
j=l 

minimizing  Jc,  one  gets  a  suboptimal  solution.  Since  this  suboptimal 
solution  depends  on  the  value  of  p,  one  would  expect  that  the  smaller  the 


value  of  p,  the  closer  the  suboptima 1  solution  should  be  to  the  optimal 
one. 

In  LQ  optimal  control,  the  intuitive  significance  of  requiring  a 
strictly  positive  input  weight  in  the  cost  function  Jc  is  to  assure  the 
boundedness  of  the  input  amplitude.  But  in  the  present  problem,  the 
input  u'(j)  equals  the  first  state  component  at  time  (j+1):  j+1)  = u(j) . 

Therefore,  the  positive  weight  q^  associated  to  ^(j+1)  would  prevent 
the  value  of  u'(j)  from  becoming  unbounded. 

Motivated  by  the  above  consideration,  we  tackle  the  problem  of 
finding  the  optimal  input  u'(n),  if  it  exists,  for  the  case  p  =  0.  In 
order  to  do  so,  we  derive  the  Riccati  matrix  for  p>0  and  observing  its 


continuous  dependence  on  p,  we  let  p-+0. 

The  Hamiltonian  has  the  expression: 

Hk  =  »s[C(k)QC(k)]+AT(k+l)[UC{k)+b1u'(k)]+|u,2(k)  (C.4) 

A(k)  =  3H/3£(k)  =  QC( k)  + UTA(k+l)  (C.5) 

0  =  3H/3u(k)  -  pu'(k)  +  b1TA(k+l)  (C.6) 

The  conjugate  equations  are: 

pC(k+l)  *  pUe(k)-b1b1TX(k+l)  (C. 7) 

A(k)  3  UTX(k+l)  +  QC(k)  (C.8) 

with  the  boundary  conditions: 

A(n)  =  Q5(n)  (C.9) 
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Now  introduce  the  matrix  P(k)  by 


A(k)  =  P(k)r(k) 

Then  equations  (C.7)  through  (C.9)  are  rewritten  as: 

[oI  +  bb1TP(k+l)]C(k+l)  =  pU£(k) 

UTP(k+lK(k+l)  =  [P(k)  -Q(k)K(k) 

P(n)  =  Q 

Upon  eliminating  £(k)  and  £(k+l)  from  the  above  equations  we  deduce: 


UTP(k+l)[I  +— -  P(k+l)]-1U  =  P(k)  -  Q 
P(n)  =  Q 


(C.10) 


Note  that 


-P12(k+1)  ....  -P1N(k+l)  ' 

P+Pu(k+1) 

0 

0 

■p+P11(k+l) 


inq  out  the  matrix  multiplication  of  (C.10)  we  deduce: 
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■  -P12(k+1)  -P13(k+1) 

P+Pl 1 ( k+l )  p+PjjCk+l) 


UTP(k+l) 


1 

0 


0 

1 


0  0 
0  0 


“P 

p+pnrk+D  0 


0 

0 


+  Q  =  P(k) 


1  0 
0  1 


(C.ll) 


The  dependence  of  the  above  Riccati  equation  on  p  is  throughout  l/[p+P11(k+l)] 
which  is  continuous  for  p>0  provided  that  P.j(k+1)^0.  That  is,  if 
P1 1  ( k )  ^  0  for  all  ke[l,n]  and  all  p>0,  we  can  let  p  +  0  in  the  above 
equation  and  deduce  the  corresponding  gain  matrix  P(k). 

Let  us  show  that  the  entry  P^(k)  never  vanishes,  i.e.,  Pjj(k)^0 
for  all  ke[l,n].  This  results  from  the  following  induction: 

p 

P(n)  =  Q  implies  that  P^(n)  =  q^ ^  =  h1  >  0  (see  Section  IV). 


From  equation  (C.ll)  we  have 

CP11(k+l)P22(k+l)  -  P122(k+1)]  +  pP22(k+l) 
pll(k)  =  qu+  p+P  1 1  ( k+l ) 

P(k+1)  being  a  semi-definite  matrix,  it  results  that: 

Pn(k+1)  >  0 
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* 


P22(k+1)  >  0 


p  / t+i \  _p  ^(k+1)  5*  0  (all  the  minors  of  the  non-negative 
HlrK  12  matrix  P(k+1)  are  non-negative.) 


Hence  it  follows  that  P^CkJ^g^^O  for  all  p>0. 

Now  we  can  let  p  =  0  in  equation  (C.ll)  and  the  resulting  equation 
(C.12)  will  generate  the  sequence  P(k+1)  corresponding  to  the  case  p-0 


i'-P12(k+l) 


-P1M(k+l) 


P(k)  =  Q  +  UTP(k+l)  0 


0  .  . 


1  .  . 


1  .  . 


2.  The  Optimum  Input  u'(j) 

The  optimum  input  u'(j)  is  derived  from  P(j+1)  by: 

u 1  ( j )  =  L(jk(j)  -  -[b1TP(J+l)b13b1TP(j+l)UC(J) 

u'(j)  =  p^l'jAT  . P1N^+1^’ 

Letting  n  go  to  infinity,  the  cost  function  J  becomes 


J  =  H  £  CT(n)Q£(n) 


Equation  (C.12),  solved  backward  for  large  n,  converges  to  a  "steady 
state"  Riccati  solution  and  the  gain  L(j)  becomes  time  invariant: 

*■  =  *-P12’  P13”"’  P1N* 

Since  only  the  first  column  of  the  matrix  P  is  of  importance  for  the  de¬ 
termination  of  C,  we  can  write  the  Riccati  recursive  equation  in  terms  of 
column  vectors  of  P. 

P  -  [pj  |PN]  and  Q  =  [qj  |qN] 


Then: 


Pj(k)  -  UTP2(kH) u'Pj(kH)  tqj 

T  Pi  l  i+1 1  t 

Pj(k>  ■  U  Pjn<ktl>-  -ifflVT)"11  Pi(ktl)  +qJ 

T  Piw(k+l)  t 

PN.i(k)  -  «V+1)-pj^U  P^lJ^H-i 
PN(k)  -  qN 
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